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SECTION  I 


INTRODUCTION 


Substantiating  postulations  through  experimentation,  or 
using  experimentation  to  draw  new  conclusions,  is 
fundamental  to  science  and  engineering.  This  service  is 
provided  to  theoretical  mechanics  by  the  field  of 
experimental  mechanics.  Before  coming  to  any  conclusions, 
experimentalists  demand  confidence  in  their  work;  as  a 
result,  many  time  consuming  experiments  are  repeated  just 
to  be  sure.  Experimental  stress  analysis  is  no  different. 
Repeating  experimental  procedures,  often  including  more  than 
one  technique,  is  the  norm  when  searching  for  the  solution 
to  engineering  problems  [1].  The  tedious  nature  of  manual 
data  extraction  and  reduction  led  researchers  to  explore 
automated  extraction  and  reduction  schemes.  Hence,  hybrid 
experimental-numerical  stress  analysis  techniques  evolved. 

The  shear  difference  method  [2]  is  an  early  example  of 
an  experimental-numerical  hybrid  method  in  stress  analysis. 
Since  then,  researchers  have  acknowledged  the  need  for 
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improved  and  more  encompassing  hybrid  techniques  [3].  Prior 
to  the  development  of  digital  image  processing,  most  hybrid 
stress  analysis  techniques  suffered  from  the  need  for  manual 
data  extraction.  Seguci  et  al.  [4]  bridged  the  gap  between 
computer  analysis  and  experimental  data  via  a  digitized 
video  picture.  Their  paper  marks  the  beginning  of  the  era 
where  researchers  started  using  digital  image  processing 
techniques  with  structural  analysis  codes  to  produce  viable 
stress  analysis  techniques  [5-6].  During  this  same  era, 
researchers  started  to  exploit  a  new  numerical  structural 
analysis  technique,  the  boundary  element  method. 

The  appearance  of  boundary  element  methods  (BEM)  was 
coincident  with  the  advent  of  digital  computers  [7],  but 
they  did  not  become  popular  until  Rizzo  [8]  introduced  the 
direct  boundary  element  method.  Still,  BEM  (sometimes 
called  boundary  integral  equation  methods)  have  yet  to 
achieve  the  acceptance  of  finite  element  methods  (FEM) .  The 
reluctance  of  the  technical  community  to  embrace  BEM  is 
largely  attributable  to  its  less  familiar  mathematical 
foundations  [9].  However,  in  the  past  ten  years,  the 
popularity  of  BEM  for  hybrid  and  purely  numerical  studies 
has  greatly  increased. 

Moslehy  and  Ranson  [10],  Umeagukwu  et  al.  [11],  and 
Balas  et  al.  [12]  introduced  early  BEM-experimental  hybrid 
techniques,  but  all  require  manual  data  extraction.  Liu 
[13]  introduced  a  hybrid  method  using  digital  image 


processing  and  elastic  BEM.  In  his  method,  digital  image 
correlation  is  used  to  find  displacements  interior  to  the 
boundary.  Liu  then  solves  the  inverse  boundary  element 
problem  to  find  the  boundary  conditions.  Liu's  work  finally 
relieved  the  experimenter  of  data  extraction 
responsibil ities . 

Hybrid  elastic  boundary  element-experimental  methods 
are  presently  an  important  research  area  [14-15].  Equally 
important  are  elastic-plastic  boundary  element-experimental 
methods,  but  there  seems  to  be  an  absence  of  research  in 
this  topic.  Elastic-plastic  boundary  element  methods 
(EPBEM)  were  first  conceived  by  Swedlow  and  Cruse  [16]  in 
1971  and  then  implemented  by  Riccardella  [17]  in  1973. 

Since  then,  EPBEM  have  grown  in  popularity.  Telles  [18] 
gives  an  excellent  historical  review  of  the  development  of 
this  numerical  method.  The  rapid  increase  in  the  popularity 
of  elastic-plastic  boundary  element  is  due  to  [18-19]: 

1.  The  BEM  solution  reduces  much  of  the  domain 
by  one  dimension.  That  is,  in  elastic 
regions,  only  the  boundary  is  discretized. 

2.  Only  the  plastic  zone  must  Le  dcr.a:-. 
discretized. 

3.  BEM  produces  a  reduced  set  of  equations 
(compared  with  FEM) . 

4.  The  user  preparation  to  run  the  codes  is 
greatly  reduced. 

5.  Expertise  requirements  are  greatly  reduced 
from  techniques  dependent  on  mesh  generation. 
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6.  Infinite  and  semi- inf inite  problems  can  be 
accurately  and  efficiently  modelled. 

7.  The  displacements,  stresses  and  strains  can 
be  selectively  calculated. 

8.  EPBEM  affords  a  high  degree  of  accuracy. 


The  EPBEM  algorithm  used  here  is  conceptually 
equivalent  to  the  rate  form  of  elastic  BEM  with  a  body  force 
[20].  The  governing  differential  equation  for  displacements 
can  be  arranged  so  that  the  homogeneous  part  is  the  Navier 
equation  for  total  displacement  rates.  The  nonhomogeneous 
part  contains  the  nonlinear  terms  associated  with  spatial 
derivatives  of  the  plastic  strain  rate.  These  nonlinear 
terms  constitute  a  pseudo  body  force.  The  boundary  element 
solution  approach  leads  to  an  equation  for  the  unknown 
boundary  conditions  in  terms  of  the  known  boundary 
conditions  (assuming  the  body  forces  are  known) .  Since  the 
plastic  strain  rates  introduce  an  additional  set  of 
unknowns,  the  boundary  element  method  requires  an  additional 
set  of  equations.  The  constitutive  model  provides  enough 
information  so  that  an  initial  strain  L2l]  solution 
procedure  can  be  successful.  This  EPBEM  algorithm  uses  the 
isotropic  hardening,  von  Mises  yield  criterion  discussed  by 
Swedlow  [22]. 
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In  the  natural  progression  of  coupling  boundary  element 
methods  an-  image  processing  techniques,  this  report 
introauces  a  two-dimensional  experimental  EPBEM  tec.in.qte. 
The  experimental  technique  used  here,  displacement  pattern 
matching,  is  a  spinoff  of  the  pattern  mapping  technique, 
first  introduced  by  Fail  [23].  Pattern  mapping  uses  digital 
image  processing  and  syntactic  pattern  recognition 
techniques  to  discern  the  motion  of  spots  applied  to  a 
specimen.  The  history  of  pattern  mapping  can  be  traced 
through  the  photogrid  method  [24-25]  and  the  fine-grid 
method  [26-27].  These  predecessors  measured  strain  by 
measuring  the  local  motion  of  large  grids  of  dots  applied  to 
a  model.  Typically,  the  measurements  were  made  by  the 
experimenter  and  thus  subject  to  a  large  degree  of  error. 
Kawasaki  et  al.  [28]  automated  the  fine  grid  method  to 
attain  good  accuracy,  but  the  equipment  is  necessarily 
complex  and  dedicated  to  the  fine  grid  task. 

•  • 

•  • 

•  • 

•  • 

•  • 

Figure  1.  Five-by-five  pattern  mapping  grid  with  an 
icon  • 
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Pattern  mapping  has  the  advantage  of  using  general 
purpose  image  processing  hardware.  Additionally,  a  regular 
grid  need  not  cover  the  entire  specimen.  In  pattern 
mapping,  one  places  an  arbitrarily  shaped  matrix  of  spots 
(for  example  a  5x5  square)  at  the  locations  in  which  the 
displacements  and  strains  are  required  (Figure  1).  Located 
in  this  matrix,  there  is  a  special  spot  (termed  an  Icon) 
that  serves  as  a  reference  spot  for  determining  the  relative 
locations  of  the  individual  spots  and  the  matrix 
orientation.  The  position  of  each  spot  in  the  matrix  is 
then  recorded  before  and  after  motion.  Differencing  yields 
the  displacement,  and  a  strain  definition  (small,  large, 
Lagrangian,  etc.)  provides  a  measure  of  strain  [23],  In  the 
proposed  hybrid  stress  analysis,  only  the  displacement  is 
required,  thus  the  name  "displacement  pattern  matching" 

(DPM) .  This  requirement  further  relaxes  the  grid  geometry 
restrictions  to  only  placing  a  spot  at  the  locations  where 
one  desires  the  displacements.  In  addition,  DPM  does  not 
use  an  Icon,  resulting  in  the  reduction  of  the  syntactic 
recognition  language  from  context  sensitive  to  context  free. 
The  combination  of  displacement  pattern  matching  for  data 
extraction  and  elastic-plastic  boundary  element  methods  for 
structural  analysis  leads  to  a  fully  automated  technique  for 
measuring  plastic  stress  and  strain. 


Three  experiments  explore  the  ability  of  the  automated 
hybrid  technique  as  a  nondestructive  stress  analysis  tool. 
All  are  plane  stress  tests  of  1100-H14  aluminum,  using  an  MTS 
testing  machine  to  apply  controlled  loads.  The  first  is  a 
simple  uniaxial  stress  test  of  a  thin  sheet.  The  second  is 
the  investigation  of  a  thin  sheet  with  a  central  hole, 
subject  to  uniform  end  loads  (perforated  strip  tension 
test);  and  the  third  is  a  thin  V-notched  sheet,  subject  to 
uniform  end  loads  (notched  tension  specimen  test) .  All 
three  experiments  serve  to  show  the  usefulness  of  this 
hybrid  DPM-EPBEM  technique  in  elastic-plastic  stress 
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SECTION  II 

THE  ELASTIC- PLASTIC  BOUNDARY  ELEMENT  METHOD 
II. 1  Linear  Elastic  BEM 


The  derivation  of  the  two-dimensional  integral 
equations  governing  linear  elasticity  and  their  solutions 
provide  an  entry  point  to  the  understanding  of  the 
underlying  concepts  of  boundary  element  methods.  The 
elasticity  equations  and  solutions  are  first  derived,  and 


then  they  are  extended  to  include  plastic  flow. 

assumptions  in  the  development  are 

The  basic 

1) 

The  problem  is  well  posed  (proper 
boundary  conditions  are  specified) . 

2) 

The  material  is  isotropic 
elastic . 

and  linear 

3) 

Small  displacement  theory 
applicable . 

is 

Elastic  boundary  element  methods  are  based  on  the 
numerical  solution  of  the  integral  equations  governing  the 
deformation  in  an  elastic  solid.  These  integral  equations 
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can  be  derived  by  starting  with  the  following  form  of  the 
virtual  work  statement: 


★ 

★ 

ti (s) Ui (s) dr+ 

t>i(P)uj.(p)dn= 

Jr  J 

n 

° i j (s,p) £ ij ( p ) dn  ; 
0 


(II-l) 


where  tj_,  Uj_,  and  bj_  denote  the  traction,  displacement  and 
body  force  vectors  respectively.  The  stress  tensor  is  a  j_  j 
and  the  strain  tensor  is  £ i j •  The  differential  dr  is  an 
element  on  the  bounding  curve  and  the  differential  dn  is  an 
interior  surface  element.  Also,  the  *-superscript  denotes 
an  arbitrary  equilibrium  set  and  no  superscript  denotes  an 
arbitrary  compatible  set  [29].  The  right  hand  side  of 
Equation  (II-l)  can  be  rewritten  using  an  identical  form  of 
virtual  work  but  with  the  significance  of  the  superscript 
and  no  superscript  interchanged.  The  resulting  equation, 
given  below,  is: 


t*(s)uj_(s)dr  + 

r 

+  bi(p)u*(p)dfi  (II-2) 

J  n 


t>i  (P)  uf  (p)  dn= 


n 


ti(s)Ui(s)dr 


This  equation  is  known  as  Betti's  Second  Reciprocal  Theorem 
[18].  In  these  two  equations,  s  represents  a  point  on  the 
bounding  curve  and  p  represents  a  point  in  the  interior. 
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Since  the  choice  is  arbitrary,  u^,  tj_,  and  bj_  are 
chosen  as  the  traction,  displacement,  and  body  force  of  the 
true  body.  The  choice  of  the  *-superscript  set  is  also 
arbitrary;  therefore,  it  is  chosen  as  Kelvin's  solution  to 
the  response  of  an  infinite,  two-dimensional,  elastic  medium 
which  is  subjected  to  arbitrarily  located,  mutually 

orthogonal  point  loads.  Kelvin's  solution  is  found  by  using 

★ 

the  Green's  function  approach  of  letting  bj_  be  the  two- 
dimensional  Dirac  delta  function  SijPj(P)  [30];  where  Pj (p) 
(for  j  =  1 , 2 )  are  two  mutually  orthogonal  unit  loads  (i=l,2). 
This  designation  is  substituted  into  the  adjoint  form  of  the 
Navier  equation  to  yield  the  following  forms  of 
displacement,  traction,  stress  and  strain: 


u*(s)  =  U*j ( s , p ) Pj (p) , 
t*(s)  =  T*j (s,p) Pj (p) , 

«ii(P)  3  Ejki(s,P) Pj (P)  ,  and 
k  i  ( P )  3  sjki(SrP) Pj (P)  > 


(H-3) 


where 


uij  <s-p>”  I^TwTg  ((3-4^)ln(r)  s  i,  -  r.ir.j, 


nj  I],!-,) G 


-  (1-2./) [rf inj-r, jni] 


»  t 


11 

Ejki  (S'P) =  i7(  1-U ) Gr  "  { (1-2l/ )  (r, k5 ij+r, j5 ik) 

-  r,  i5  j k+2r .  ir,  jr/k)  , 

( I I -4 ) 

sjki  (S'P)=  4^  ( i-u  )  r  ^  1-2./ )  (ri  k5  i  j+r ,  j  5  ki-r ,  i5  j  k) 

*  2r, ir , j r , k ) ; 

In  Equations  (II-4),  r2  =  (Xp-Xs)2  +(Yp-Ys)2,  i/  is  Poisson's 
ratio,  n^  is  the  outward  nomal  vector  to  the  bounding 
surface,  dr/dn  is  the  normal  derivative,  and  G  is  the  shear 
modulus.  The  function  U*ij(s,p)  is  known  as  Kelvin's 
fundamental  solution;  where  T*ij(s,p),  E*jki(s»P)»  and 
S*jki(s,p)  are  the  traction,  strain,  and  stress  respectively 
that  are  associated  with  the  fundamental  solution. 

it 

Sokolnikoff  [31]  provides  a  detailed  development  of  U^j  and 

* 

most  boundary  element  texts  include  a  derivation  of  T j_  j  , 

★  ★ 

Ejki<  anc*  sjki  (f°r  example  see  Ref.  [20]). 

Somigliano's  Identity  for  displacements  [32]  is  the 
basis  for  the  boundary  element  method.  This  identity  is 
found  by  substituting  the  fundamental  solution  into  Betti's 
Second  Reciprocal  Theorem. 


ui(p)= 


Uij ( s , p ) t j (s)dr- 


Tij ( s , p) u j (s)dr 


+ 


Uij (s) bj (p)dfi 


n 


( 1 1  -  5 ) 
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Equation  (II-5)  can  be  restated  as,  given  all  the  boundary 
conditions  (the  body  force  is  always  considered  known) , 
information  at  any  point  inside  the  boundary  is  attainable 
without  finding  the  entire  domain  solution.  In  other  words, 
selective  calculation  of  interior  displacements  is  possible 
by  a  line  integral  around  the  bounding  curve.  An 
appropriate  real  estate  analogy  is  that  a  buyer  need  never 
see  the  property,  he  only  needs  to  walk  the  perimeter  to 
know  every  detail  of  the  lot  and  the  house.  This  is  a  truly 
remarkable  statement.  The  solution  to  Navier's  equation  is 
reduced  in  dimension  (in  the  sense  that  only  the  boundary 
must  be  examined) .  The  reduced  dimensionality  (surface  to 
line  integral)  makes  BEM  a  very  attractive  numerical 
solution  technique.  But  still,  the  unknown  boundary 
conditions  remain  unknown,  so  the  solution  is  not  yet 
tractable . 

The  above  dilemma  is  easily  remedied  by  taking  the 
limit  of  Equation  (II-5)  as  the  field  point,  p,  approaches 
the  surface  [8].  Mathematically  stated  as 

Lim  Uj_(p),  where  t  is  a  point  on  r.  (II-6) 

p-  i 

Since  somewhere  on  the  boundary  r=0  (i=s),  and  since 
T*ij(s,p)  is  proportional  to  1/r  and  U*ij  is  proportional  to 
ln(r),  the  limit  in  Equation  (II-6)  produces  a  singular 
integral.  This  singular  integral  exists  [8]  and  is 
interpreted  in  the  principal  value  sense  [33].  The 


existence  proof  of  the  limit  in  Equation  (II-6)  follows  by 
segmenting  the  integration  path  to  isolate  the  singular 
point,  then  taking  the  limit  of  the  integral  as  the  singular 
integration  path  shrinks  to  zero. 

The  existence  of  the  right  hand  side's  second  integral 
in  Equation  (II-5)  is  verified  by  the  approach  described 
above.  Using  the  notation  in  Figure  2, 


■  r 

*  * 

lim  Ti j (s , i ) uj ( s ) dr  =  lim  Tjj  (s, l ) Uj (s) dr 

p-*o  Jp  p-oJr 

(H-7) 

,  .  * 

+  lim  Tji (s, 2 ) ujGdr ' . 

p  -*  0  p  _  p  ' 


The  first  integral  on  the  right  hand  side  of 
Equation  (II-7)  exists  in  the  ordinary  sense  but  the 
integral  around  the  singular  circuit  must  be  viewed 
carefully.  The  singular  integral's  existence  is  confirmed 
by  making  the  following  two  substitutions: 

1)  T*ij  =g (i)/p:  where  g (?)  is  order  one  and  2)  dr=pds 


’«2  f*2 

lim  g(0)d$  =  g [8)  d 9.  (II-8) 

p-  0  S  9 

1  1 

Generally,  the  limit  in  equation  (II-7)  is  not 
evaluated  until  assumptions  about  the  interpolation  of  Uj 
over  the  boundary  r  are  made  (see  section  III.  4).  Lachat 
[34],  following  closely  the  arguments  of 
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Figure  2.  Circuit  around  the  singularity. 


Muskhelishvili  [33],  showed  that  Equation  (II-7)  evaluates 
to  C^jUj  regardless  of  the  interpolation  function.  In  the 
strict  sense  Cj_j  is  a  constant  and  is  dependent  on  surface 
smoothness;  but  in  the  numerical  solution,  C^j  is  also 
dependent  on  the  interpolation  function. 

The  same  arguments  hold  for  the  other  two  integrals  on 
the  right  hand  side  of  Equation  ( 1 1  —  5 ) .  In  the  first 
integral  U*ij  is  proportional  to  ln(p)  and  in  the  third 
integral  the  area  element  is  p&pdt  .  In  both  cases,  the 
limit  evaluates  to  zero.  Assembling  the  aforementioned 
results  yields  the  boundary  element  constraint  equation. 


Cij (s)uj (s) 


Ti.j  (s,  t )  uj  (s)dr  = 

r 


Uij  ( s ,  2  )  t j  (s)  dr 

r 


+ 


U*j ( s , p) b j (p)dn. 

n 


( II-9 ) 
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For  a  well  posed  problem,  Equation  (II-9)  provides  the 
unknown  boundary  conditions.  These  boundary  conditions,  m 
turn,  can  be  used  in  Equation  (II-5)  to  give  the 
displacement  at  any  desired  interior  point.  As  an  example, 
if  the  surface  traction  is  prescribed  on  the  surface,  then 
Equation  (II-9)  yields  the  unknown  surface  displacements. 
Then  Equation  (II-5)  gives  the  displacement  solution  in  the 
domain.  This  procedure  is  valid  for  prescribed  tractions, 
prescribed  displacements,  or  mixed  boundary  conditions. 

II. 2  Elastic-Plastic  BEM 

The  governing  equations  in  elastic-plastic  boundary 
element  methods  are  developed  along  the  same  lines  as 
elastic  BEM;  the  main  difference  is  the  idea  of  initial 
strain.  This  section  will  describe  the  initial  strain 
concept  and  its  effect  on  the  BEM  development. 

Additionally,  the  need  for  a  second  equation  is  illustrated 
and  a  plastic  flow  rule  is  introduced  to  meet  that  need. 

Swedlow  and  Cruse  [16]  were  first  to  formulate  the 
initial  strain  boundary  element  theory.  Since  then, 
Riccardella  [17] ,  Mendelson  and  Albers  [35],  Mukherjee  [36], 
Telles  [18]  and  others  have  developed  initial  strain 
algorithms  to  solve  engineering  problems.  Although  initial 
strain  theory  is  not  restricted  to  a  particular  flow  rule, 
it  will  be  described  in  the  context  of  isotropic  work 


hardening  and  incremental  plasticity.  In  accordance  with 
accepted  terminology  and  notation,  rate  and  increment  will 
be  used  interchangeably  and  superior  dots  denote  time 
differentiation. 

In  an  elastic-plastic  body,  strains  can  be  categorized 
as  either  elastic  or  plastic.  In  initial  strain  EPBEM ,  the 
plastic  strain  is  assumed  to  act  as  a  residual  strain;  each 
lead  increment  is  just  a  deviation  from  that  residual 
strain.  The  adjective  "residual"  refers  to  the  idea  that 
the  solid  knows,  by  some  internal  mechanism,  the  plastic 
strain  increment  before  the  application  of  the  load 
increment.  When  the  increment  is  applied,  the  material  only 
deforms  elastically.  Loosely  interpreted,  the  initial 
strain  concept  states  that  plastic  deformation  occurs  before 
elastic  deformation.  This  idea  of  initial  strain  is 
analogous  to  Martin's  discussion  of  convergence  by  the 
superposition  of  elastic  and  residual  stresses  [37].  As  an 
aside,  initial  strain  has  little  meaning  in  finite  strain 
plasticity  since  numerically,  this  is  a  load  and  unload 
process.  Initial  strain,  however,  is  acceptable  in 
incremental  plasticity. 

The  derivation  of  the  integral  equations  governing 
elastic-plastic  BEM  follow  the  same  procedure  outlined  in 
Section  II. 1.  This  time,  however,  the  starting  point  is  a 
rate  form  of  virtual  work: 
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★  • 

r  *  .  r 

tj_  (s)  u  ( s)  dr  + 

fc>i  (P)  Ui  (p)  dn  = 

r 

o 

<7*j  (s,p)  £  ij  (P)dn  ; 
a 


( 1 1  - 1  o ; 


The  superior  dot  denotes  time  differentiation,  the  terms 
designated  by  the  *-superscript  denote  an  arbitrary 
equilibrium  set,  and  the  terms  with  no  superscript  denote  an 
arbitrary  compatible  set.  The  total  strain  tensor  can  be 
decomposed  into  elastic  and  plastic  parts;  thereby, 
rewriting  Equation  (11-10)  as 


f  *  • 

★ 

ti(s)ui(s)dr+ 

fc>i(p)ui(p)dn= 

r 

n  J 

*  *  P 

° i j (s,p) £ ij (p)dfi 

n 


★  .  e 

*ij (s,p) £ ij (p)dn. 

JO 


(ii-ii) 


During  each  load  increment  both  elastic  and  plastic 
deformation  are  present.  Even  though  plastic  deformation 
occurs,  the  stress  rate  is  still  related  to  the  elastic 
strain  rate  by  Hooke's  law  [38].  Accordingly,  the  stress 
rate  can  be  found  by  isolating  the  elastic  part  of  the  total 
strain  increment  and  using  Hooke's  Law: 

*ij  -  2G(£ij  -  «?j)  +  ^dkk  -^k)*  (H-12) 


Equation  (11-12)  allows  the  following  substitution  to  be 

*  •  e  *  • 

made  in  equation  (11-11)  :  ^ij^ij  =  «ij®ij* 

Equation  (11-11)  can  then  be  written  in  a  form  similar  to 
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Equation  (II-2)  by  applying  the  rate  form  of  virtual  work 

(with  the  significance  of  the  superscripts  interchanged] 

★  • 

to  the  e  1  j cr  i  j  term; 


*  • 

• 

*  ★ 

tj_(s)Uj_(s)dr  - 

tf  (s)  Uj_  (s)  dr  + 

r 

r 

( P )  u  i  ( p )  da 


(H-13) 


r.  f  *  -p 

-  bi(p)ui(p)dn  =  crfj  (S,p)  f  ij  (p)dfi. 

.  a  J  a 

The  difference  between  this  integral  equation  and 
Equation  (II-2)  is  the  integral  which  includes  the  plastic 
strain  rate  term.  As  in  Section  II. 1,  the  terms  with  the 
♦-superscript  are  chosen  as  Kelvin's  solution. 

(Equation  ( 1 1  —  3 ) )  and  the  terms  with  no  superscript  are 
chosen  as  the  displacement  rate,  traction  rate,  body  force 
rate  and  strain  rate  of  the  true  body.  These  designations 
lead  to  Somligliano1 s  Identity  for  inelastic  materials; 


Uj  (s) 


U*j (s,p) ti(s)dr 

r 


u*j  (s,p)t>i(p)dn  + 

r 


( s , p ) u^(s)dr 


4- 


f 

Ja 


(11-14) 

*  .P 

o ij (s,p) € ij (p) da  . 


The  limiting  process  to  find  the  elastic-plastic  boundary 
element  constraint  equation  is  identical  to  that  described 
in  Section  II. 1.  The  plastic  strain  rate  integral  exists 
and  has  no  contribution  to  j  . 
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Cij(s)  Uj(s)  +  T*j  (s,i)ui(i)dr  =  f  u*j(s,i)  tj_  ( ^  )  dr 

J  n  J  q 

+  u*j  (s,^)bi(i)dn  +  [  s*ki(s,  i)  e^fijdn; 

.  r  Jr 

where  tf*ki(s,p)  =  S* j ki (s , p) Pj (p)  is  used. 

If  the  plastic  strain  is  known  apriori  and  if  the 
problem  is  well  posed,  then  Equation  (11-15)  yields  the 
unknown  boundary  conditions.  The  boundary  information  could 
then  be  used  to  find  the  displacements  in  an  elastic-plastic 
solid  (Equation  (11-14)).  The  problem  here  is  that  prior 
knowledge  of  the  plastic  strain  rate  is  just  an  idealization 
of  the  initial  strain  solution  approach.  In  reality,  the 
plastic  strain  rate  is  an  unknown  quantity;  therefore,  an 
additional  relation  is  needed. 


;iastic-Plastic  Flew  Rule 


An  isotropic  work  hardening,  von  Mises  flow  rule  is 
introduced  to  satisfy  the  above  requirement  for  an 
additional  relation.  Although  a  multitude  of  flow  rules 
have  been  proposed  [39],  an  isotropic  work  hardening  flow 
rule  provides  an  easily  understood  constitutive  model.  It 
also  provides  a  reasonable  approximation  to  material 
response  while  loading  takes  place.  However,  one  should  be 
aware  that  the  von  Mises  yield  criterion  is  limited  in  that 
it  does  not  provide  for  hysteresis  [4uj. 


Figure  3.  Isotropic  hardening  yield  surface. 

The  yield  function  of  an  isotropic  hardening  material 
mathematically  describes  a  yield  surface  that  maintains  its 
shape  but  increases  in  size  as  a  function  of  the  plastic 
deformation.  As  an  example  (Figure  3),  f(g)  defines  the 
yield  surface,  where  g  is  some  measure  of  plastic 
deformation.  The  loading  path  AB  is  purely  elastic,  where  B 
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is  an  infinitesimal  amount  inside  the  yield  surface.  As  the 
solid  is  loaded  along  BC,  plastic  deformation  commences  and 
the  yield  surface  expands.  Now,  unloading  along  CD  is 
purely  elastic  until  D  is  reached.  As  a  result,  the  elastic 
domain  is  enlarged. 

There  are  infinitely  many  forms  of  isotropic  work 
hardening  to  choose  from;  Koiter  [41]  presents  a  general 
development  of  work  hardening.  Hill  [42]  presents  a  detailed 
discussion  of  work  hardening  and  its  applications,  and 
Malvern  [43]  gives  a  concise  review  of  the  prevalent  models. 
The  isotropic  work  hardening  rule  used  here  employs  a  von 
Mises  yield  criterion;  where  the  yield  function  is  single 
valued  and  independent  of  hydrostatic  pressure.  The  yield 
criterion  may  be  stated  as 


f  (<7ij)-F(J2)>0; 


(11-16) 


where  f  is  the  yield  function  and  J2  is  the  second  invariant 
of  the  stress  deviatoi. 

Stability  in  Drucker's  sense  leads  to  the  statement  of 
normality  [44]  which,  in  turn,  leads  to  a  functional  form  of 
the  proposed  constitutive  model.  Normality  is  expressed  as 


•p.  _  'LL 
^  -  A  a  a 


(H-12) 
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where  a'  is  a  proportionality  factor  that  depends  on 
J2  and  dJ2  [45].  Once  on  the  yield  surface, 
f  ( cr  i  j  )  =  F  ( J2 )  <  additional  plastic  deformation  requires 
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f (<* ij ) >F ( J2 ) •  Therefore  one  may  write 

‘P  . d f  a f  •  .  , 

fij  -  A.  “mn1  (11-18) 

5a i j  3c7ran 

where  a  is  a  different  proportionality  factor  and  is  a 
function  of  J2  [45].  Normality  produces  a  linear  relation 
between  plastic  strain  rate  and  the  stress  rate.  In  order 
to  make  Equation  (11-18)  useful  for  computations,  specific 
assumptions  need  to  be  made  about  the  form  of  F(J2).  Under 
the  assumption  of  von  Mises  yielding,  the  yield  function  is 
the  von  Mises  equivalent  stress  (creq)  and  is  defined  as  [43] 

aeq  =  -/3  Sij  Sij  ;  (11-19) 

where  S]_j  =  -  5  ijo)c]<  is  the  stress  deviator.  The  3/2 

factor  in  Equation  (11-19)  ensures  ae q  equals  the  yield 
strength  when  uniaxial  load  is  applied.  For  a  work 
hardening,  von  Mises  material,  an  equivalent  plastic  strain 
rate  can  be  defined; 


=  J 2  «P. 

T 


(11-20) 


such  that  dWp  =  aeq<eq  where  dWp  is  the  plastic  work 
increment  [40].  The  multiplication  of  Equation  (11-18)  by 
a ij  and  substitution  of  dWp  =  ae q«eq  [40]  into  the  result, 
leads  to  the  following  expression  for  A. 


A 


(H-21) 
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d  f  2 

Then,  by  noting  —  damn  =  dJ2  and  3J2  =  ae f  the 

dc,mn 

proportionality  factor  is  simplified  to 


4a 2  Em 
eq  1 


(11-22 ; 


where  Ex  is  the  tangent  modulus  [46]  of  the  equivalent 
stress-equivalent  strain  curve.  The  definitions  of 
equivalent  stress  and  strain  are  such  that  under  uniaxial 
loading  ET  is  also  the  tangent  modulus  of  the  uniaxial 
stress-strain  curve. 

Finally,  substituting  Equation  (11-22)  into  (11-18) 
leads  to  a  specific  form  of  the  isotropic  flow  rule,  which 
gives  a  relation  between  the  plastic  strain  rate,  stress 
deviator  and  the  stress  rate. 


9  sij  smn  -a 

4  E«-  a  2  1 

^  eq 


(11-23) 


Elastic  deformation  is  included  by  combining 

Equation  (11-23),  Hooke's  law,  and  the  total  strain  rate 
e  p 

tensor  (d« ^ j  +  d<ij)  to  yield  a  relation  for  the  total 
strain  rate 


?ij  ~  2  g  ^kk  1  + 


9  Sii  Smna 

4  Em  o 2 
1  eq 


mn; (11-24) 
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and  the  stress  rate 


° ij  -  2G[t  i j 


l-2i 


')  5  i  j  ‘  Ick  ]  ~ 


3G  Sii  5 


ran  f  mn 


° gq  ( 1+Et/ 3G) 


■;  (11-25) 


where  incompressibility  ( s' ^ k=0  )  i-s  used.  By  substituting 
Equation  (11-25)  into  Equation  (11-23)  the  final  form  of  the 
von  Mises,  isotropic,  work  hardening  constitutive  model  is 
obtained . 


•P 

4  ij 


3  sij  smn  ‘mn 
2  a  2  (1+  ET/3G) 


(11-26) 


II. 4  Numerical  Implementation  of  EPBEM 

Simple  domains,  simple  boundary  conditions  and 
knowledge  of  the  loading  history  are  all  needed  if  exact  • 
solutions  to  the  elastic-plastic  boundary  element  equations 
(Equations  11-14,  11-15  and  11-26)  are  to  be  found.  Some 
exact  solutions  are  documented  [47],  but  seldom  are 
researchers  so  fortunate. 

This  section  will  discuss  an  iterative  solution  to 
the  elastic-plastic  boundary  element  problems.  In  doing  so, 
the  surface  and  domain  discretization,  the  interpolation 
functions  used  to  approximate  the  traction  rate  (tjj. 
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► 


displacement  (UjJ  rate,  and  the  plastic  strain  rate  (,'P^j) 
are  presented.  The  system  matrix  formulation  is  discussed 
in  detail.  This  discussion  includes  the  calculation  of 
diagonal  elements,  the  calculation  of  surface  stresses,  and 
the  handling  of  geometric  discontinuities.  Finally,  the 
iterative  solution  scheme  for  the  elastic-plastic  problem  is 
detailed.  In  order  to  simplify  the  following  discussion, 
body  forces  are  neglected.  Also,  for  convenient 
referencing.  Equations  (11-14),  (11-15),  and  (11-26)  are 

rewritten  below. 


Uj(p)  = 


Uij (s,p)ti(s)dr 


Tij (s,p)ui(s)dr  + 


*  .P 

Sjki(s'P) f  ki (P) dn ■ 


Q 


(11-27 .A) 


Cij(s)  u i  ( s )  + 


Tj_ j  (sri)Qi(i)dr  = 


Uij(s,i)  ti (i ) dr 


s jki (s>  ^ )  *  ki  ( ^  ^  * 


(11-27. B) 


3Sii  Smn 
2  U+  ET/3G) 


I 


(11-27. C) 
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The  daily  advances  made  in  memory  and  speed  of  modern 
computers  amaze  even  the  most  avid  computer  enthusiasts. 
However,  today's  computers  have  nowhere  near  the  capacity  to 
exactly  model  the  functions  and  domains  commonly  encountered 
in  many  mechanics  problems.  For  this  reason,  domain 
discretization  and  function  interpolation  have  become  a 
science  unto  themselves.  Boundary  element  methods,  just 
like  finite  element  methods  (FEM)  and  finite  difference 
methods  ( FDM) ,  need  a  certain  amount  of  surface  and  domain 
quantization.  Fortunately,  the  integral  approach 
significantly  reduces  the  amount  of  discretization 
necessary.  As  seen  in  Equations  (11-27) ,  both  surface  and 
domain  integrals  must  be  evaluated.  But  in  BEM,  the  domain 
integrals  are  only  evaluated  in  the  plastic  zone  and 
therefore  represent  a  significant  advantage  over  FEM  and 
FDM.  Another  important  difference  is  that  the  domain 
discretization  is  used  only  to  evaluate  the  domain  integral. 

A  further  difference  between  finite  element  methods  (finite 
difference)  and  boundary  element  methods  lies  in  the 
solution  approach.  FEM  and  FDM  approximate  the  solution 
over  the  domain,  so  the  solution  is  highly  dependent  on  the 
quantization  and  function  interpolation.  In  contrast,  the 
boundary  element  governing  equation  is  solved  exactly  and 
all  approximations  are  made  in  the  boundary  conditions. 


2' 


k 


Figure  4.  Boundary  modelled  with  straight  line 
segments . 

Some  cases  have  been  reported  where  BEM  has  afforded  as 
much  as  an  order  of  magnitude  improvement  in  accuracy 
over  FEM  [48]  . 

The  boundary  element  method  requires  some  way  of 
describing  the  bounding  curve  of  a  region;  here,  boundary 
discretization  is  achieved  with  straight  line  segments 
(Figure  4).  Straight  line  segments  can  simply  and 
accurately  approximate  a  boundary;  in  many  cases  (plates, 
beams,  etc.)  the  approximation  is  exact.  If  necessary,  the 
surface  can  be  more  closely  modeled  by  increasing  the  number 
of  segments.  It  should  be  noted  that  there  is  a  practical 
limit  to  the  size  and  number  of  these  straight  line 
segments,  but  that  limit  is  rarely  encountered. 
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Now  assumptions  must  be  made  about  the  behavior  of  the 
traction  and  displacements  over  these  straight  line 
segments.  By  definition,  the  tractions  are  derived  from  the 
spatial  derivatives  of  displacement.  Although  consistency 
in  derivatives  is  often  ignored,  Brebbia  and  Walker  [49] 
postulated  and  Liu  et  al.  [6]  later  showed  that  a  consistent 
boundary  element  solution  greatly  exceeds 


Figure  5.  Representative  boundary  element. 

comparable  inconsistent  routines  in  both  computational  speed 
and  accuracy.  In  addition,  when  only  displacement  boundary 
conditions  are  prescribed,  upwards  to  a  50%  increase  in 
speed  can  be  achieved  [6].  This  last  attribute  is  of 
particular  interest  in  hybrid  experimental-numerical 
methods . 

The  proposed  EPBEM  algorithm  strives  for  a  balance 
between  simplicity  and  accuracy.  To  this  purpose,  the 
boundary  displacements  are  assumed  to  vary  linearly  along 
the  boundary  segments.  Correspondingly,  the  tractions  are 
assumed  constant  over  each  segment.  In  BEM  jargon  these  are 
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called  constant  (traction)  and  linear  (displacement) 
elements.  Specifically,  the  boundary  is  broken  up  into  N 
elements,  each  element  is  joined  to  its  neighboring  elements 
at  nodes.  The  linear  variations  between  nodes  can  be 
described  by  (Figure  5) 

ui(x)=Ui(l+2s/L)+Ui(l-2s/L)  for  -L/2  <  s  <  L/2;  (11-28) 

where  the  superscripts  1  and  2  denote  node  1  and  2 
respectively,  and  the  subscript  i  is  1  for  x-displacement 
and  2  for  y-displacement .  The  form  of  Equation  (11-28) 
dictates  that  the  displacement  information  is  located  at  the 
nodes,  whereas  the  constant  element  traction  information  is 
located  at  the  mid  points  of  each  segment.  As  an  aside, 
more  sophisticated  discretization  schemes  are  available 
[13,54] . 

The  plastic  strain  rate  still  must  be  quantized 
before  the  numerical  procedure  can  be  described.  The 
plastic  zone  is  discretized  much  like  FEM.  Triangular  cells 
are  chosen  with  the  plastic  strains  assumed  to  be  constant 
over  each  cell.  This  may  seem  identical  to  FEM  but  in 
contrast  to  FEM,  the  cells  are  used  only  to  evaluate  the 
domain  integral.  If  the  interpolating  functions  are  given, 
then  the  constraint  equation  (II-27.B)  and  interior 
displacement  equation  (11-27. A)  can  be  rewritten  as 
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(11-29. B) 


Equation  (11-29. A)  represents  a  system  of  2N  equations 
for  2N  unknowns.  Its  implementation  is  best  described 
through  Figure  6.  Each  set  of  two  rows  produced  by  Equation 
(11-29. B)  is  a  result  of  placing  the  source  point  (xs,  y5) 
at  the  center  of  an  element  and  then  integrating  its 
influence  around  the  surface.  The  source  point  is  theji 
moved  from  element  to  element  until  the  circuit  is  complete. 
Two  rows  are  produced  because  Equation  ( 11-29. B)  represents 
both  x  and  y  directions.  The  above  procedure  leads  to  a 
system  of  equations  in  the  form  of 


[T](u>=[U](t)+[S] (<P} ; 


(11-30) 
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Figure  6.  Depiction  of  BEM  solution  methodology. 

where  [T]  is  the  displacement  coefficient  matrix,  [U]  is  the 
traction  coefficient  matrix,  (u)  are  the  boundary 
displacement  rates  (some  of  which  are  known) ,  and  (t)  are 
the  boundary  traction  rates  (some  of  which  are  also  known) . 
The  plastic  strain  coefficient  matrix  is  [S],  and  («P) 
contains  the  plastic  strain  rates  (all  of  which  are 
unknown) .  Grouping  the  known  boundary  conditions  and 
plastic  strain  rate  information  on  the  right  hand  side 
reduces  Equation  (11-30)  to 

(A] (v}=(q)+[S] (<p) ;  (II-31) 
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where  [A]  is  the  system  matrix,  {v}  contains  the  unknown 
boundary  conditions,  and  (q)  contains  the  known  boundary 
conditions  multiplied  by  the  appropriate  coefficients.  As 
will  be  discussed  later.  System  (11-31)  can  now  be  solved 
for  the  unknown  boundary  conditions. 

There  is  little  difficulty  in  using  a  Gaussian 
numerical  quadrature  [50]  to  evaluate  most  of  the  line 
integrals.  However,  care  must  be  exercised  when  evaluating 
the  elements  containing  the  singularity  (diagonal  elements) 
as  well  as  evaluating  the  domain  integrals.  The  diagonal 
elements  exist  but  are  singular?  therefore,  either 
analytical  or  special  numerical  techniques  are  necessary. 
The  domain  integral  takes  special  care  since  there  is  a 
certain  amount  of  controversy  over  the  derivative  of 
Equation  ( 11-29. B). 

The  first  alternative  in  evaluating  the  singular 
diagonal  elements  is  to  surround  the  singular  point  with 
many  small  sub-elements  and  then  integrate  numerically  [51] 
The  second,  and  preferred,  method  is  to  integrate 
analytically  along  each  singular  segment  [17,35].  In  any 
case,  the  principal  value  coefficient  matrix  (C^j)  is 
calculated  analytically  [17]  (shown  below)  and  added  to  the 
diagonal  terms  to  form  the  system  matrix. 
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The  domain  integration  scheme  used  in  this  algorithm  is 
chosen  to  avoid  the  current  controversy  over  the  correct 
expression  for  the  spatial  derivative  of  Equation  (11-27. A). 
This  derivative  is  a  dominating  aspect  since  it  defines  the 
strain  within  the  body.  The  ongoing  debate  revolves  around 
the  existence  and  the  proper  definition  of  the  derivative  of 
the  plastic  strain  rate  integral.  Although  well  documented 
by  Telles  [18]  in  his  introductory  chapter,  the  salient 
features  of  the  controversy  are  presented. 

The  strain  at  any  point  in  the  domain  is  defined  by 
2t  =  (ui,j  +  uj,i)-  Many  early  researchers  failed  to 
recognize  that  the  integration  path  around  the  domain 
singularity  changes  with  the  load;  this  path  change 
prevents  pulling  the  derivative  directly  under  the  plastic 
strain  rate  integral  [52-55].  Bui  [56]  was  the  first  to 
recognize  this  error  and  formulated  the  correct  derivative; 
Telles  and  Brebbia  [57]  expanded  Bui's  ideas  and  then  later 
carried  them  out  [21].  The  above  techniques  are  both 
elegant  and  independent  of  the  plastic  strain  rate  shape 
function  but  they  are  sometimes  difficult  to  use. 
Specifically,  high  aspect  ratio  cells  cause  the  integration 
scheme  to  give  erroneous  results.  Since  the  EPBEM  algorithm 
presented  here  is  of  the  constant  cell  type,  the  general 
approach  championed  by  Telles  and  co-workers  will  be 
foregone  in  lieu  of  an  implementatior.ally  easier  technique. 
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Figure  3.  Scheie  used  in  domain  integration. 

The  assumption  that  plastic  strain  rate  is  constant 
over  each  interior  cell  is  not  only  computational ly  simple, 
it  also  makes  the  above  discussion  inconsequential.  The 
constant  cell  approximation  enables  direct  evaluation  of  the 
plastic  strain  rate  integral  [17,36,37]  by  subdividing  the 
cell  (Figure  8).  The  integrals  over  triangle  1  and  triangle 
2  are  algebraically  added  to  the  integral  over  triangle  3  to 
give  the  desired  integral  over  the  cell  (triangle  4). 
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The  unit  directions  are  defined  so  that  the  algebraic  signs 
are  correct.  Since  the  plastic  strain  rate  kernal  is 
singular  at  the  base  point,  a  small  region  surrounding  this 
point  is  excluded  in  the  kernal  evaluation  (Appendix) .  This 
integration  scheme  is  equally  valid  when  the  source  point 
coincides  with  the  load  point  (Figure  9).  No  great 
difficulty  exists  in  evaluating  the  integrals;  and  the 
question  surrounding  the  existence  of  the  derivatives  is 
moot  since  plastic  strain  rates  are  found  via  direct 
differentiation . 
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Figure  9.  Singularity  inside  the  integration 
cell . 

Another  point  of  computational  interest  is  the 
evaluation  of  stress  rates  along  the  boundary.  Since  Sj^i 
has  a  singularity  0(r-2),  the  integral  approach  to  calculate 
stress  rates  near  the  boundary  introduces  large  errors 
(since  r<<!  near  the  boundary).  Instead,  the  EPBEM  employs 
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the  technique  sketched  by  Hartmann  [58].  Although 
introduced  in  the  context  of  elastostat ics ,  the  method  is 
equally  valid  for  elasto-plasticity .  There  are  seven 
unknown  quantities  on  any  smooth  boundary  point  in  a  two- 
dimensional  body,  three  stress  rates  j  and  four 
derivatives  u^j.  There  are  also  seven  equations  available: 

1)  Hooke's  Law 

'a  ij  =  2Ge  ij  +  ~  T  ^  i  j  fkk  (11-37. A) 

1  —  Zv 

2)  The  definition  of  traction  rate 

tij  =  a  j  nj  (11-37. B) 


3)  The  derivative  of  u^  along  the  arclength. 


(11-37. C) 


Note  that  Equation  (11-37. A)  uses  only  the  elastic  portion 
of  the  strain  rates  so  as  to  keep  the  number  of  unknowns  at 
seven;  Equation  (11-37. c)  uses  only  the  elastic 
displacements.  Since  displacements  are  approximated  by 
linear  functions,  the  derivatives  along  the  arc  length 
(Equation  (11-37. C)  are  found  directly.  When  constant 
elements  approximate  displacements,  finite  difference 
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methods  are  typically  used  to  provide  the  arclength 
derivatives  [59].  The  seven  equations  assembled  in  matrix 
form  appear  as 


1 

0 

0 

a 

0 

0 

-X 

*11 

0 

nl 

0 

"2 

0 

0 

0 

0 

*  2  2 

^1 

0 

n2 

nl 

0 

0 

0 

0 

*  11 

^2 

0 

1 

0 

—  X 

0 

0 

a 

uei,l 

= 

0 

0 

0 

0 

-n2 

0 

nl 

o 

^S2,  1 

•  p 

u  l,s 

0 

0 

1 

0 

-G 

-G 

0 

•  p 

U®1,2 

0 

0 

0 

0 

0 

-n2 

0 

"1  . 

#  p 

L  u  2,2 

.  u02,s  . 

where  the  traction  rates  and  the  elastic  tangential 
displacement  rate  derivatives  are  a  product  of  the  boundary 
element  constraint  equation,  x  ■  2Gis/{l-2v),  a  -  -x-2v  ,  nj_ 
are  the  normals  and  G  is  the  shear  modulus.  The  equations 
for  stress  rates  are  calculated  with  Cramer's  rule  and 
listed  below: 
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where 

A  =  x  2n^  +  n]_-a2nxn2  ^  , 

Bi  =  an]_2  +  \ri22,  B2  =  nltl-n2t2' 

Ci  =  A2-a2,  C2  =  u2 ' snl~ul , sn2 , 

D  =  Gan^+(2GA  -a2-A  2  )  r»x2n|^  Gan^  , 

E  =  Ari]_2+an2  . 

With  the  preliminary  numerical  eccentricities  taken 
care  of,  it  is  now  time  to  present  the  final  part  of  the 
iterative  EPBEM  solution  puzzle.  First  the  boundary 
integral  equation  for  the  total  strain  rate  is  needed.  By 
substituting  Equation  (11-27. A)  into  strain-displacement 
equation,  the  total  strain  rate  looks  generically  like 


(11-40) 
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sijra  (s,p)dfl  is  given  m  the  Appendix. 

As  previously  mentioned,  the  singular  nature  of  Sjki(s-P) 
prevents  bringing  the  derivative  directly  under  the 
integral;  therefore  Sjk;j_(s,p)  is  analytically  integrated  and 
then  differentiated  (Appendix) . 

The  iterative  elasto-plastic  boundary  element  solution 
technique  is  similar  to  the  successive  approximation  scheme 
used  by  Roberts  and  Mendelson  for  stress  function  solutions 
[60].  The  iterative  approach  used  here  is  less  cumbersome 
since  the  total  strain  rates  are  found  directly.  The 
instructions  below  list  the  steps  in  the  EPBEM  solution 
technique  and  Figure  10  provides  a  flow  chart  for  quick 
referencing. 
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1)  Apply  an  elastic  load  to  bring  the  most 
highly  stressed  cell  to  a  stress  state  near 
yield. 

2)  Find  the  total  elastic  stress  and  strain  at 
each  cell.  For  each  cell,  use  the  elastic 
strain  as  the  initial  guess  for  accumulated 

and  use  elastic  stress  as  the  initial 
guess  for  the  current  total  stress. 

3)  Apply  a  load  increment  large  enough  to  cause 
incipient  yielding.  Use  the  load  increment 
and  the  initial  guess  for  the  plastic  strain 
rate  in  the  following  form  of  Equation  (II- 
31)  to  find  the  unknown  boundary  conditions. 

[A] { v}=  ( q}  +  6( q)  +  [AZ](tp> 

6 (q)  -  load  increment 

4)  Use  (v)  from  step  (3)  to  find  the  total 
strain  rates  at  each  cell  with 
Equation  (11-40) . 

5)  Calculate  the  elastic  strain  rate  at  each 
cell  with  e  ei j  —  € ^ j  —  €  Pi ^ . 

6)  Use  Hooke's  Law  (11-12)  to  calculate  the 
stress  rate  at  each  cell. 

7)  Add  the  stress  increment  to  the  current  total 
stress  in  each  cell  and  use  it  to  find  the 
corresponding  tangent  modulus  (Et)  from  the 
uniaxial  stress  strain  curve. 

8)  Check  each  cell  to  see  if  it  is  beyond  yield; 

&  e  f  >  ?  ¥  * 

where  cref  is  the  von  Mises  equivalent  stress 
(also  the  yield  function)  and  Y  is  the 
current  uniaxial  yield  strength  for  the  cell. 

9Y)  If  the  cell  has  yielded;  then  calculate  the 
new  plastic  strain  rate  of  the  cell  with  the 
isotropic  work  hardening,  von  Mises  flow  rule 
(11-33. C).  Use  the  new  total  stress,  Et  from 
step  (7)  and  the  current  plastic  strain  rate 
in  Equation  (11-33. C). 
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9N)  If  the  cell  is  still  elastic,  set  <Pj_j  equal 
to  zero. 

10)  Check  for  eP^j  convergence  against  the 
previous  value. 

1 IN)  If  the  tolerance  is  not  met,  subtract  the 
stress  rate  found  in  step  (6)  from  the 
current  total  stress  found  in  step  (7) .  Then 
go  to  step  (5)  and  start  again. 

11Y)  If  the  error  tolerance  is  met,  approximate 
the  new  uniaxial  yield  stress  for  each 
yielded  cell  by  the  current  von  Mises 
equivalent  stress. 

12)  If  the  problem  is  done,  quit.  If  the  problem 
is  not  finished  go  to  step  (3)  and  start  a 
new  increment  cycle. 


With  possible  exception  of  step  (11N),  the  instructions 
listed  above  are  self-explanatory.  For  step  (11N),  the 
current  guess  for  the  stress  rate  must  be  subtracted  from 
the  known  previous  total  stress  after  each  pass  through  an 
unsuccessful  iteration.  It  is  improper  to  add  each  current 
guess  for  the  stress  rate  to  the  last  increments  total 
stress  plus  the  last  iterations  stress  rate  guess.  The 
current  total  stress  would  then  be  the  sum  of  the  previous 
total  stress  and  all  the  stress  rate  guesses. 

The  above  iterative  solution  to  elasto-plastic  boundary 
element  problems  is  reasonably  efficient  and  has  acceptable 
accuracy.  The  technique  is  sufficiently  general  so  that  it 
is  easily  extended  to  other  than  boundary  element  methods. 
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Figure  10. 


Flow  Chart  of  the  iterative  elasto-plast ic 
boundary  element  solution  technique. 
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When  large  scale  yielding  occurs,  there  is  a  possibility 
that  the  plastic  zone  will  exceed  the  celled  region.  if 
such  a  case  arises,  one  expects  a  reduction  in  EPBEM 
accuracy . 
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SECTION  III 

DISPLACEMENT  PATTERN  MATCHING 

The  spread  of  image  processing  applications  in 
photomechanics  is  a  consequence  of  the  availability  of 
equipment  at  affordable  costs.  Hardware  prices  continue  to 
decline  and  the  power  of  personal  computers  continues  to 
increase  Economical  mass  storage,  networking,  graphics, 
and  application  specific  software  all  have  increased  the 
proliferation  of  image  processing  in  the  photomechanics 
field  [61].  For  the  most  part,  image  processing  te^.niques 
in  photomechanics  are  grouped  into  two  categories,  fringe 
analysis  [62-65]  and  correlation  techniques  [6,66-68].  A 
third  category,  pattern  mapping  [23,69],  is  emerging  as 
alternative  with  many  attractive  qualities.  Pattern  mapping 
is  a  non-destructive  technique  measuring  displacement  and 
strain  and  is  based  on  image  processing  and  syntactic 
pattern  recognition.  It  is  a  process  by  which  "black"  spots 
on  a  "white"  background  are  followed  in  a  double  exposure 
type  scheme.  This  technique  offers 
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1)  Full  field  measurement . 

2)  Sub-pixel  registration. 

3)  Greatly  reduced  CPU  time. 

4)  A  highly  automated  technique. 

5)  Relative  ease  of  specimen  preparation. 

Pattern  mapping  was  originally  developed  as  a  general 
purpose  strain  measurement  technique  capable  of  discerning 
large  rigid  body  rotation  and  translations.  It  also 
included  a  variety  of  strain  definitions.  Fail  [23] 
suggested  that  reductions  in  run  time  were  achievable  by 
tailoring  the  pattern  mapping  to  a  specific  class  of 
problems.  Displacement  pattern  matching  ( DPM)  is  a  result 
of  this  suggestion.  The  displacement  pattern  matching 
algorithm  is  developed  by  assuming  that  the  experimenter  can 
prevent  appreciable  rigid  body  motion  and  that  displacements 
are  the  only  information  of  interest.  As  a  result,  a  highly 
automated,  very  fast  displacement  tracking  scheme  emerged. 

The  DPM  algorithm  can  be  broken  into  three  blocks: 
image  segmentation,  feature  extraction,  and  matching. 

These  three  terms  are  common  pattern  recognition  jargon  but 
may  be  new  to  the  experimental  mechanics  community.  To 
clarify  the  transition  into  the  new  terminology,  formal 
definitions  of  image  segmentation,  feature  extraction, 
matching  and  a  fourth  term,  registration,  are  provided. 


47 


Image  Segmentation  is  the  separation  of  the  image 
into  different  regions,  each  having  specific  and 
related  properties.  It  is  the  first  step  in  all 
pattern  recognition  schemes  and  attempts  to 
describe  or  classify  the  image  [70, 71],  An 
example  might  be  to  describe  a  landscape;  that  is, 
segmentation  of  a  satellite  photograph,  where 
classification  might  entail  distinguishing  rivers 
from  the  background  [72]. 

Feature  Extraction  is  the  quantification  of 
certain  characteristics  of  an  image  which  are 
usable  in  further  image  information  classification 
[73].  This  process  discerns  important  foreground 
features  from  the  background.  Examples  might  be 
lengths,  areas,  edges  or  curves  [74]. 

Matching  refers  to  a  class  of  operations  of 
comparing  to  two  or  more  images  [75].  Matching  is 
typically  used  in  time-varying  imagery,  where 
motion  detection  is  simply  a  matter  of  detecting 
differences  in  successive  images  [76,77]. 

Registration  is  the  comparison  of  two  images  of  a 
scene  taken  from  different  perspectives. 
Registration  is  used  as  a  measure  of  resolution 
attainable  with  matching  algorithms  [78]. 


Capturing  a  digital  image  and  the  equipment  required  to  do 
so  are  briefly  discussed  in  the  next  section.  The  ensuing 
sections  describe  the  segmentation,  feature  extraction,  and 
matching  blocks  of  the  DPM  algorithm.  The  final  section 
verifies  DPM  through  rigid  body  motion  tests  and  will 
attempt  to  assess  the  important  characteristics  that  control 
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The  human  eye  scans  the  analog  intensity  distribution 
which  is  created  by  light  reflecting  from  a  scene  and  sends 
this  information  to  the  brain  to  be  interpreted.  Upon 
receiving  the  signal,  the  brain  processes  the  scene  and 
classifies  the  image.  Ideally,  computer  vision  should  be 
capable  of  emulating  this  process.  A  quantization  of  the 
intensity  distribution  is  necessary  if  a  computer  is  to  have 
any  hope  of  mimicking  the  recognition  processes  of  the  human 
brain.  The  most  primitive  medium  of  communication  for  a 
computer  is  numbers  (digits).  Therefore,  it  is  logical  that 
an  efficient  quantization  technique  would  be  able  to  convert 
the  analog  intensity  distribution  to  a  digital 
approximation,  hence  the  term  digitize. 

A  digital  image  is  a  numerical  approximation  of  a  real 
scene,  with  local  intensity  distributions  approximated  by 
gray  levels  [79].  Gray  levels  are  integers  ranging  from 
zero  to  an  upper  limit  set  by  the  hardware.  The  eight  bit 
digitizing  board  used  in  these  experiments  affords  0-255 
gray  levels.  These  gray  levels  represent  the  darkness  and 
are  ordered  to  form  a  matrix  representation  of  the  image. 


An  example  of  the  imaging  procedure  is  given  in  Figure  11 
where  a  black  spot  and  its  digital  representation  are 
provided.  Key  traits  to  note  include:  black  appears  as 
lower  gray  levels  than  white,  the  two  pixel  (picture 
element)  fuzzy  transition  region  between  black  and  white 
[80],  and  the  noise  in  the  image.  When  looking  at  scenes 
containing  high  contrast  sub-regions,  the  human  eye  has 
little  difficulty  locating  boundaries  between  light  and 
dark.  For  computer  vision  hardware,  this  process  is  not  so 
simple.  One  would  expect  sharp  contrast  borders  in  the 
scene  to  be  represented  by  a  jump  in  gray  level,  when  in 
fact,  the  imaging  system  sees  a  finite  transition  between 
light  and  dark.  This  fuzzy  transition  is  a  product  of 
hardware  limitations  in  the  camera  and  in  the  digitizing 
board.  The  camera  sensing  array  is  divided  into  a  finite 
number  of  sensing  cells.  Cell  crosstalk,  frequency 
limitation  characteristics,  and  cell  overlap  are  typical 
contributions  of  the  camera  to  fuzzy  regions  [81-82],  In 
addition,  the  mismatch  between  the  camera  array  size  and  the 
digitizing  hardware  array  size  increases  the  fuzziness.  In 
general  these  array  sizes  differ  considerably.  A  typical 
example  would  be  four  digitizing  pixels  representing  three 
camera  pixels,  i.e.  a  smearing  effect.  Although  the 
mismatch  is  not  usually  this  dramatic,  its  effect  is 


noticeable . 
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Figure  11.  A  black  spot  and  its  digital 
representation . 
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In  Figure  11  the  spot  is  uniformly  black;  unfortunately 
noise  degradation  prevents  the  image  from  reflecting  the 
uniformity.  Degradation  is  apparent  in  the  fluctuating  gray 
level  representation  of  black  in  the  spot  and  results  from 
optical  or  electronic  sources.  Spherical  aberration,  coma, 
astigmatism,  non-uniform  illumination,  and  dust  are  just  a 
few  examples  of  optical  noise  [S3].  Electronic  noise  may 
come  from  thermal  effects  in  the  imaging  array,  the 
digitizing  process,  electro-magnetic  interference,  or  may 
just  be  inherent  to  the  imaging  array  [80-82].  Regardless 
of  the  source,  DPM  regards  noise  as  random  and  therefore 
having  a  zero  mean  distribution. 

Now  that  some  of  the  ideas  and  important  aspects  of  a 
digital  images  have  been  illustrated,  the  equipment  and 
procedure  for  obtaining  a  digital  image  can  be  described. 
Figure  12  depicts  a  typical  personal  computer  based 
digitizing  system.  The  imaging  device  used  here  is  a  Sony 
XC-38  CCD  camera,  a  Datacube  IVG128  Video  Acquisition  and 
Display  Board  (384H  x  512V)  which  quantizes  the  image  with 
its  digital  form  displayed  on  a  PVM-1271Q/1371QM  Sony  studio 
monitor.  All  image  processing  is  done  on  a  PC's  Limited 
286-8  Personal  computer  equipped  with  a  math  co-processor 
and  EGA  graphics. 
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Figure  12.  Typical  personal  computer  based  digitizing 
system . 


Digitizing  is  software  controlled.  On  command,  the 
frame  buffer  stores  the  image  as  a  matrix  of  signed  bytes 
representing  the  gray  levels.  This  two-dimensional  matrix 
is  also  available  to  the  user  on  software  control.  On 
software  control,  the  analog  signal  from  the  CCD  (charge 
couple  device)  camera  is  quantized  with  an  analog-to-digital 
converter  and  then  mapped  to  the  appropriated  gray  levels 
by  a  software  selectable  look-up  table.  The  image  size 
stored  in  the  frame  buffer  is  384H  x  485V  to  conform  with 
the  RS-170A,  60Hz  United  States  interlaced  video  standard. 
The  384H  x  512V  image  conforms  to  the  CCIR,  50Hz.  European 
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interlaced  video  standard-  The  frame  buffer  is  accessible 
to  the  CPU  (central  processing  unit)  through  the  personal 
computers  I/O  page  [84].  VNA  Systems  Incorporated's  Digital 
Correlation  Metrology  software  package  provides  software 
control  of  the  digitizing  process. 


ill .  2  Segmentation 


A  decided  advantage  of  displacement  pattern  matching  is 
the  limited  prior  knowledge  required  for  success.  Success 
of  the  technique  hinges  on  the  placement  of  high  contrast 
spots  (dark  spots  on  a  light  background)  at  locations  where 
displacement  is  desired.  Implicit  in  this  knowledge  is  the 
idea  that  the  images  are  binary,  that  is,  the  gray  levels 
should  be  either  bJack  or  white  [85].  In  keeping  with  the 
automated  methodology,  DPM  does  not  need  prior  knowledge 
about  the  number  of  spots,  spot  area  or  any  other 
distinguishing  feature.  The  primary  function  of 
segmentation  is  to  locate  the  important  regions  of  an  image. 
Since  DPM  images  are  binary,  segmentation  consists  of 
identifying  those  individual  regions  defined  as  black  (the 
spots)  since  they  contain  the  desired  information  [86]. 

Segmentation  typically  (as  is  done  here)  starts  with 
some  form  of  characteristic  feature  thresholding  [87].  In 
essence,  a  value  of  some  characteristic  feature  of  the  image 
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is  chosen  as  a  basis  on  which  to  make  decisions.  This 
pivotal  value  is  known  as  the  "threshold".  In  DPM,  the 
characteristic  feature  is  the  gray  level  and  the  threshold 
is  the  parameter  used  to  discriminate  between  black  and 
white.  All  gray  levels  above  the  threshold  are  white  and 
all  gray  levels  below  are  black.  The  transformation  to  a 
true  binary  gray  scale  enables  a  border  following  scheme  to 
locate  the  spots.  Discerning  if  each  spot  is  a  true  spot  or 
if  it  should  be  discarded  is  segmentation's  culminating 
decision. 

Once  the  threshold  is  chosen,  raster  scanning  starts. 

A  raster  scan  is  nothing  more  than  a  row  by  row  or  column  by 
column  evaluation  of  each  pixel  for  some  prescribed 
property;  in  this  case  to  find  the  first  spot  (gray  level  < 
I-j>  ),  When  DPM  locates  a  black  pixel,  the  gray  levels  of 
its  eight  nearest  neighbors  are  checked.  If  at  least  one 
neighbor  is  black  then  DPM  has  located  a  spot.  If  none  of 
its  neighbors  are  black,  then  DPM  considers  the  center  pixel 
as  a  false  spot.  Once  a  true  spot  is  located,  a  general 
border  following  scheme  developed  by  Rosenfeld  [88]  and 
specialized  to  only  find  outside  edges  [89,23]  'frames'  the 
spot.  Framing  the  spot  is  the  process  by  which  a  complete 
clockwise  circuit  around  the  border  is  traversed,  storing 
only  the  left  most,  right  most,  highest,  and  lowest  pixel 
locations  (L,  R,  T,  B  respectively)  for  the  feature 
extraction.  To  account  for  the  transition  between  black  and 
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white,  the  frame  automatically  expands  by  two  pixels  in  all 
directions  (Figure  11) .  With  the  first  spot  framed,  the 
raster  scanning  continues.  This  process  is  repeated  until 
the  entire  image  is  examined. 

False  spots  larger  than  one  pixel  in  diameter  are 
entirely  possible,  especially  when  the  scene  is 
significantly  magnified.  An  automated  method  of  discerning 
false  from  true  spots  is  therefore  required.  Experience 
indicates  that  the  areas  of  false  spots  are  much  smaller 
than  their  legitimate  cousins  (area  is  defined  as 
A  -  [ R-L] [ B-T] ) .  A  spot  area  larger  than  one  third  the 
average  spot  area  (averaged  over  every  framed  spot)  is  a 
true  spot.  This  definition  is  completely  arbitrary  but  it 
proved  to  be  99%  reliable  for  the  spot  sizes  used  in  this 
research . 


III. 3  Feature  Extraction 

Features  are  selected  measurements  which  are  either 
invariant  or  less  sensitive  to  typical  system  distortions 
[93];  they  tend  to  distinguish  the  object  from  its 
background.  Two  types  of  features  are  present  in  the  DPM 
images,  statistical  and  structural  [90].  The  statistical 
features  are  the  frame  boundaries:  L,  R,  T,  B,  and  the  spot 
centroid.  The  structural  features  are  the  spots  themselves. 
Since  segmentation  simultaneously  extracts  the  frame  and  the 
spots,  only  the  centroids  need  be  extracted. 


The  spot  centroids  are  calculated  using  the  true  gray 
levels  in  the  frame  and  pixel  locations  (i,j)  in  the  matrix 
[80].  Mathematically  stated  the  centroids  are 
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where  I(i,j)  is  the  gray  level  at  element  (i,j)  in  the  image 
matrix.  The  influence  of  the  zero  mean  noise  on  the 
centroid  calculations  is  presumed  to  average  to  zero. 


III.  4  Image  Matchin* 


Zero  rigid  body  motion  reduces  the  syntactic  pattern 
recognition  language  [91]  of  DPM  from  context  sensitive  to 
context  free  [92].  Pattern  mapping  requires  each  spot  to  be 
located  at  the  same  "address"  before  and  after  deformation 
[23],  In  computational  space,  this  means  each  spot  must 
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always  be  surrounded  by  the  same  neighbors.  If  this  is  not 
the  case,  that  spot  is  out  of  context  which  means  /[t  makes 
no  syntactic  sense.  DPM  makes  no  such  requirement.  The 
spot  locations  are  stored  in  the  order  in  which  they  are 
found.  Pattern  matching  [93]  is  then  used  to  find  the  spots 
in  all  ensuing  registered  images.  Ordering  the  spots  by 
their  correct  "addresses"  is  no  longer  crucial  to  finding 
the  next  image's  spots;  therefore,  the  spots  do  not  need  to 
be  in  context  to  make  syntactical  sense. 

An  instructive  approach  to  describing  pattern  matching 
is  to  take  a  figurative  step  back  and  simply  look  at  the 
physical  nature  of  a  typical  solid  mechanics  problem.  The 
question  is  "what  are  the  displacements  at  points  a,  b,  and 
c?"  To  answer  this  question  correctly,  each  spot  must  be 
examined  after  deformation  and  compared  to  its  previous 
location.  Successive  load  increments  produce  a  sequence  of 
comparisons.  In  elasticity  and  incremental  plasticity, 
assuming  spot  motion  to  be  small  relative  to  spot  size  is 
restrictive  but  certainly  reasonable.  This  indicates  that 
each  sequential  image  should  closely  resemble  its 
predecessor.  The  above  description  of  events  is  not  unique; 
Rosenfeld  and  Kax  [88]  (pages  50-51)  give  an  almost 
identical  description  of  events  but  their  discussion 
pertains  to  time-varying  imagery. 
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Time-varying  imagery  is  a  popular  form  of  pattern 
matching  [75] .  It  is  a  sequential  registration  technique 
used  in  motion  detection,  object  tracking,  and  dynamic  scene 
analysis;  therefore,  it  is  well  suited  for  discerning 
displacements.  Here,  spot  motion  is  measured  by  local 
matching  of  the  segmented  spots.  In  essence,  the  location 
of  a  spot  in  the  image  is  assumed  near  its  location  in  the 
previous  sequential  image.  Therefore,  instead  of  raster 
scanning  the  entire  new  image,  only  the  previous  framed  area 
is  scanned. 

If  the  displacement  field  carries  the  spot  out  of  the 
original  frame  area  then  the  search  region  is  expanded  to 

An  = [ L-R ] [ B-T ] n2 ,  n=2 , 3 , 4 . . .  (III-4) 

where  An  is  the  area  of  the  search  region  and  n  is  the 
number  of  search  region  expansions  required  to  find  the  spot 
(Figure  13).  When  the  search  region  encompasses 
more  than  one  spot  (Figure  14)  or  a  wrong  spot  only 
(Figure  15),  a  one-to-one  match  is  not  guaranteed  and  the 
matching  process  fails.  Therefore,  over-expansion  defines  a 
restriction  on  DPM,  but  at  the  same  time  somewhat  relaxes 
initial  assumptions  of  the  matching  scheme.  The  spot 
displacement  must  be  smaller  than  or  equal  the  distance  from 
the  original  spot  to  the  closest  spot  in  the  current  image 
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Figure  13.  Search  region  expansion. 

(Figure  15) .  This  critical  distance  is  easily  gaged.  For 
example,  in  a  square  array  of  spots,  the  maximum 
displacement  can  never  be  greater  than  0.7071  times  the 
smallest  original  spot  spacing.  Spot  displacements  are  no 
longer  restricted  to  being  small  compared  with  the  spot 


diameter. 


spot  a  w 
ower  image 


spot. 


60 


3PUT  IN  THE 
PREVIOUS  IMAGE 
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Figure  15.  Search  region  encompassing  the  wrong  spot. 
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A  special  search  region  expansion  case  occurs  near  the 
image  edge.  When  the  displacement  field  carries  a  spot 
close  to  the  perimeter  of  the  image,  the  expansion  may  try 
to  extend  the  search  region  beyond  the  image  border,  the 
results  of  which  are  unpredictable.  The  search  region  is 
therefore  contoured  to  the  image  edge  if  the  region 
encounters  the  perimeter. 

With  a  one-to-one  correspondence  established  between 
segments  in  successive  images,  statistical  feature 
extraction  can  proceed.  The  centroid  of  each  spot  is 
calculated  by  Equations  (III-2)  and  (III-3)  and  the  x-  and 
y-coordinates  of  the  centroids  in  sequential  images  are 
subtracted  to  find  the  displacements 

ux  -  Xn-Xn_!  (III-5A) 

and 

uy  =  Yn-Yn-i-  (III-5B) 

In  these  displacement  definitions,  n  refers  to  the  current 
image.  The  units  of  displacement  are  pixels  but  can  be 
converted  to  a  units  of  length  through  calibration. 

III. 5  Rigid  Body  Motion  Tests 

In  this  closing  section,  the  limitations  of 
displacement  pattern  m  itching  are  explored.  A  series  of 
rigid  body  motion  tests,  each  with  different  spot  diameters, 
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evaluate  the  accuracy  of  DPM.  In  the  process,  general 
accuracy  characterist ics ,  accuracy  versus  spot  size,  effect 
of  image  averaging  on  accuracy  and  the  consequences  of 
search  region  over-expansion  are  examined. 

Each  rigid  body  motion  test  consists  of  tracking  the 
5x5  array  of  spots  shown  in  Figure  13  through  four 
horizontal  0.025  (±0.001)  inch  increments.  Spot  areas  are 
easily  altered  by  magnification  with  a  zoom  lens. 

These  tests  produce  results  that  are  difficult  to 
assimilate  and  make  system  evaluation  unnecessarily  hard.  A 
statistical  approach  reduces  the  information  to  a  more 
understandable  form.  For  a  given  magnification,  the  25 
spots  in  the  array  are  subjected  to  four  equal  in'-^ments , 
yielding  100  displacement  values  which  should  be  identical. 
Unfortunately,  the  displacements  returned  by  DPM  are  not  all 
equal.  For  the  sake  of  comparison,  the  standard  deviation 
of  the  100  displacement  values  is  calculated  and  normalized 
by  the  accompanying  mean.  This  number  is  a  measure  of  the 
amount  of  deviation  from  the  applied  constant  displacement 
and  serves  as  a  vehicle  to  evaluate  spot  size. 

The  results  in  Table  1  show  that  there  is  a  definite  lower 
limit  on  spot  size.  Although  difficult  to  say  precisely,  a 
spot  diameter  smaller  that  ten  pixels  is  considered  suspect. 
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Table  1.  DPM  error  as  a  function  of  spot  size. 


DIAMETER 

(PIXELS) 

STANDARD  DEVIATION 
DIVIDED  BY  THE  MEAN 

6 

0.019872 

9 

0.003052 

I  • 

0 . 003203 

16 

0 . 004553 

18 

0.004842 

24 

0 . 007174 

Ideally,  all  spots  in  the  array  should  give  the  sane 
displacement  (for  rigid  body  motion).  System  and 
experimental  error  degrade  this  idealization  and  produce  a 
certain  amount  of  variance  from  spot  to  spot.  The  effect  of 
system  error  can  be  seen  by  examining  one  increment  of  the 
16  pixel  diameter  spots.  The  average  displacement  (over  the 
25  spots)  is  5.183  pixels  with  a  standard  deviation  of 
0.0039  pixels.  Averaging  images  from  the  same  scene  reduces 
the  influence  cf  random  system  noise  (cell  crosstalk, 
thermal,  electromagnetic,  etc.).  To  show  this,  four  images 
of  each  16  pixel  diameter  spot  scene  are  averaged  before 
using  DPM.  The  results  show  a  significantly  reduced 
standard  deviation  (0.00255  pixels)  which  suggests  a  more 
accurate  measure  cf  displacement  (5.181  pixels). 

Henceforth,  all  hybrid  experiments  use  four  i.mage  averages. 

"False  system  displacement"  is  a  result  of  the  variance 
in  displacement  between  spots  in  the  same  image  and  can  be 
evaluated  through  sel f- vey istration .  If  an  image  is 
registered  with  itself,  each  spot  should  measure  zero  rigid 
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body  motion.  But  instead  the  system  displacements  are 
returned.  Table  2  reveals  the  displacement  errors 
introduced  by  the  system  and  represents  the  lower  accuracy 
limit  afforded  by  DPM  and  the  hardware. 


Table  2.  DPM  system  displacement  as  a  function  of 
spot  size. 


DIAMETER 

(PIXELS) 

MEAN  SYSTEM 
DISPLACEMENT 
(PIXELS) 

6 

0 . 002251 

9 

0 . 001435 

14 

0. 001063 

16 

0.001632 

18 

0.002117 

24 

0.001562 

The  above  errors  are  inherent  to  feature  extraction  and 
are  somewhat  abstract;  that  is,  no  definite  values  can  be 
used  to  discriminate  good  from  bad.  Search  region  over¬ 
expansion  is  different  in  this  respect;  over-expansion 
defines  the  limit  to  the  displacement  increment  size.  For 
horizontal  rigid  body  motion,  one  expects  the  matching  to 
fail  when  the  displacement  reaches  one  half  the  spot 
spacing.  This  limit  is  easily  verifiable  through  a  simple 
experiment.  A  5x5  array  of  0.05-inch  diameter  spots  is 
again  tracked  through  0 . 025-(±0 . 001)  inch  increments;  but 
this  time  the  increments  are  continued  until  the  matching 
fails.  The  results  are  shown  in  Table  3  for  a  spot  spacing 
of  0.2  inches  and  a  spot  diameter  in  pixels  of  14.  As 
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predicted,  the  matching  is  precise  until  a  displacement 
greater  than  0.1  inch  is  applied.  The  limit  on  displacement 
increment  will  not  be  as  easily  identifiable  when  structural 
distortions  are  encountered,  but  a  safe  rule  of  thumb  is  not 
to  exceed  one  fifth  of  the  smallest  spot  spacing. 

I 

Table  3.  DPM  search  region  over  expansion. 


APPLIED 

MEASURED 

DISPLACEMENT 

DISPLACEMENT 

( inches) 

( inches) 

0.0 

0.00005 

0 . 025 

0.0249 

0.05 

0 . 0499 

0 . 075 

0 . 0751 

0 . 1 

0.0996 

0.125 

0.0286 

Displacement  pattern  matching  is  an  accurate,  efficient 
means  of  tracking  displacements.  All  DFM  calculations  are 
performed  on  an  80287-8  based  personal  computer  with  run 
times  averaging  about  31  seconds  per  increment. 
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SECTION  IV 

HYBRID  DPM-EPBEM  TECHNIQUE 


Both  the  elastic-plastic  boundary  element  and 
displacement  pattern  matching  algorithms  have  been 
individually  verified,  but  the  ultimate  goal  here  is  to 
construct  a  single  successful  stress  analysis  tool  b\ 
meshing  these  two  techniques.  Accordingly,  this  chapter 
presents  a  hybrid  DPM-EPBEM  and  its  experimental 
verification.  The  following  sections  describe  the 
experimentation,  specimen  preparation,  and  verification 
tests . 


IV. 1  Working  Environment 

In  general,  combining  two  programs  written  in  the  same 
language  is  a  trivial  matter.  But,  when  the  two  programs 


are  written  on  two  different  computers  and  use  two  differen 
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Figure  16.  Graphical  Working  Environment 

languages  the  process  is  significantly  more  complicated. 
Here,  DPM  is  written  on  a  PC's  Limited  286-8  personal 
computer  and  EPBEM  is  written  on  a  Digital  VAX  11/750. 

These  two  computers  have  no  direct  means  of  communicating 
with  each  other,  so  a  third  program,  PROCOMM,  is  used  as  an 
interpreter.  To  run  DPM-EPBEM,  the  image  of  each  increment 


is  digitized  and  the  resulting  digital  representation  is 
stored  in  such  a  way  that  DPM  can  either  access  these  files 


68 


interactively  or  in  batch  mode.  In  any  case,  the  user 
controls  the  program  through  a  graphical  working 
environment . 

The  working  environment  shown  in  Figure  16  provides 
critical  information  for  real  time  evaluation  of  the 
program's  progress.  A  binary  pseudo-image  shows  the  scene 
currently  being  analyzed.  Included  in  this  picture  is  a 
window  encompassing  the  segment.able  portion  of  the  image. 
The  pseudo-image  is  produced  by  averaging  each  set  of 
sixteen  nearest  neighbors  and  then  thresholding.  The  resul 
is  a  94H  x  121V  matrix  representation  of  a  binary  image. 
This  representation,  as  opposed  to  a  384H  x  484V,  is 
necessitated  by  the  limited  dynamic  memory  afforded  by 
Microsoft's  DOS.  On  the  right  of  the  pseudo-image  is  a 
histogram  of  the  true  image  with  the  location  of  the 
threshold  indicated  by  the  vertical  line.  The  correspond  in 
value  of  the  threshold  is  given  below  the  histogram.  Also 
provided  below  the  histogram  is  the  percentage  of  pixels 
with  the  most  frequent  gray  level.  All  this  information 
aids  in  evaluating  the  performance  of  the  DPM . 

Above  the  pseudo  image  and  histogram  shown  in 
Figure  16,  is  a  region  of  the  screen  reserved  for  any 
necessary  software  commands.  That  portion  of  the  screen  is 
called  the  editor  zone.  After  each  set  of  displacements  is 
determined  by  DPM,  the  communications  program,  PROCOMM,  is 
evoked  through  a  batch  file.  The  batch  file  logs  onto  the 
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VAX,  transfers  the  displacement  files,  and  starts  EPBEM. 
During  this  process,  the  editor  zone  displays  each  batch 
command  as  it  is  executed.  Once  the  VAX  attains  control, 
all  of  the  graphics  are  temporarily  suppressed;  and  after 
EPBEM  has  completed  the  necessary  calculations ,  the  VAX 
maintains  control.  VAX  control  enables  the  user  to  transfe 
the  stress  and  strain  information  back  to  the  host  computer 
(personal  computer) ,  to  return  to  the  host  computer,  or  to 
log  off.  Regardless  of  which  is  chosen,  all  of  the  desired 
stress  information  is  stored  on  the  VAX  and  can  easily  be 
retrieved  for  further  analysis. 


IV. 2  Experimental  Setup  and  Specimen  Preparation 


This  section  provides  a  description  of  the  equipment 
and  of  the  general  experimental  setup.  Since  EPBEM  is  a 
two-dimensional  code,  only  plane  stress  or  plane  strain 
loading  conditions  can  be  examined.  All  specimens  arc  of 
the  plane  stress  type  and  are  produced  from  a  ductile 
aluminum  alloy  (1100-H14).  Loads  are  applied  by  a  Material 
Testing  System  (MTS)  machine.  A  schematic  of  the  setup  is 
given  in  Figure  17.  The  testing  machine  is  manually 
controlled  to  produce  either  constant  load  or  constant 


stroke . 


i  MTS  C:NTHCL-53 
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Figure  17.  Schematic  of  the  experimental  set  up 
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As  described  in  Section  III,  DPM  depends  on  the 
application  of  high  contrast  spots  to  the  specimen.  The 
aluminum  specimens  are  painted  flat  white  in  order  to 
produce  the  high  contrast  background  for  the  black  spots. 
Initial  painting  attempts  resulted  in  the  paint 
chipping  and  cracking  under  lead.  Since  DPM  perceives 
chips  or  cracks  as  spots,  one  must  undertake  means  tc  step 
this  phenomenon  from  occurring.  Each  specimen  is  first 
thoroughly  sand  blasted  to  remove  gross  deposits  of  aluminum 
oxide.  They  are  then  cleansed  in  a  bath  of  isopropyl 
alcohol.  Following  air  drying,  each  specimen  is  scrubbed 
carefully  with  AMCHEM's  Alumiprep  33  (stock  no.  DX533)  to 
remove  the  remaining  corrosive  oxidation  and  to  chemically 
etch  the  surface.  The  final  treatment  is  a  generous  scrub 
and  bath  in  AMCHEM's  Alodine  1201  (stock  no.  DX503).  This 
treatment  chemically  stabilizes  the  aluminum  surface  and 
provides  a  long  lasting  paint  adhesive.  An  additional 
advantage  of  this  procedure  is  that  it  leaves  a  distinctive 
amber  stain  where  the  coating  process  is  successful.  The 
specimens  are  then  air  dried  and  painted  with  flat  white 
KRYLON  enamel.  The  entire  process,  from  start  to  finish, 
requires  about  ten  minutes  per  specimen  and  provides  a 
painted  surface  which  does  not  chip  or  crack  until  near 
catastrophic  failure  loads  are  applied. 

Spot  application  followed  an  iterative  process  until 
finding  a  satisfactory  procedure.  The  primary  concerns  here 
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are  uniformity  in  spot  diameter  and  darkness.  Best  results 
are  obtained  by  using  typical  "rub-ons"  commonly  found  in 
office  supply  stores;  the  grammatical  periods  are  far  more 
consistent  in  area  and  darkness  than  any  attempted  painting 
or  staining  procedures.  The  periods  are  simply  rubbed  cn  at 
the  locations  where  the  displacement  is  desired. 

Three  plane  stress  experiments  are  presented  in  order 
of  increasingly  complex  states  of  stress.  The  first  test; 
examined  is  a  thin  sheet  under  uniaxial  tension,  the  second 
test  is  a  thin  rectangular  sheet  with  a  central  hole, 
subjected  to  tensile  end  loads  (perforated  strip  tensile 
test)  and  the  final  test  is  a  thin  rectangular  sheet  with 
symmetric  ninety  degree  V-notches.  The  I100-H14  aluminum 
used  for  all  specimens  has  the  following  material 
properties:  E=10,713.  kpsi,  ^=.33,  and  Yo=15,0Q0  psi. 

EPBEM  requires  the  uniaxial  stress  strain  curve  in  the 
iterative  solution.  To  accommodate  EPBEM,  the  uniaxial  true 
stress-engineering  strain  curve  is  found  experimentally  by 
applying  known  loads  to  a  uniaxial  test  specimen  with  the 
MTS  machine  and  monitoring  a  one-inch  gage  length  with  the 
MTS  supplied  extensometer .  The  stress-strain  data  is  fitted 
with  two  piecewise  continuous  curves.  The  elastic  and  knee 
portion  are  fitted  with  the  Ramberg-Osgood  Law  [94];  and  the 
work  hardening  portion  is  approximated  with  a  linear  fit. 

The  approximate  and  experimental  stress-strain  curves  are 
given  in  Figure  18. 
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Figure  13.  Uniaxial  stress-strain  curve 
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Prior  to  each  experiment,  the  digitizing  system  is 
calibrated  carefully  so  that  DPM  results  are  accurately 
converted  from  pixels  to  units  of  length.  The  calibration 
constant  is  a  function  of  magnification  so  that  the  first 
step  in  calibrating  the  system  is  to  set  the  desired 
magnification  and  to  keep  it  fixed  for  ail  load  increments. 
Next,  a  5x5  array  of  black  dots  is  placed  in  the  focal  plane 
and  given  five  equal  vertical  displacements;  the  array  is 


/  -t 


digitized  after  each  step.  Since  all  of  the  displacements 
are  equal,  a  linear  relation  exists  between  the 
displacements  returned  by  DPM  (in  pixels)  and  the  measured 
displacements;  the  slope  of  this  curve  is  the  calibration 
constant.  The  steps  to  calculating  the  calibration  constant 
are  as  follows: 

1)  Apply  each  load  increment. 

2)  Store  the  four  image  averages  of  each 
scene . 

3)  Calculate  the  average  rigid  body 
displacement  of  the  25  spots  in  each 
scene . 

4)  Plot  the  average  displacement  (in  pixels) 
versus  the  measured  displacement  (in  units 
of  length)  . 

5)  Use  a  linear  least  squares  fit  to  find  the 
slope . 

The  imaging  array  of  the  digitizing  board  used  in  these 
experiments  is  rectangular;  therefore,  the  y-cal ibration 
constant  differs  from  the  x-cal ibration  constant.  To 
accommodate  the  rectangular  imaging  array,  the  above 
calibration  procedure  is  repeated  with  horizontal 
displacements.  A  typical  experiment  produced  x-  and  y- 
calibration  constants  of  369.5  and  523.3  pixels/inch, 
respectively.  with  system  resolution  on  the  order  of 
.1  pixels,  the  displacement  resolution  is  on  the  order  of 


.0001  inch. 
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I V . 4  Uniaxial  Tension  Test 


Before  describing  the  hybrid  experiments,  it  should  b 
noted  that  no  attempt  is  made  to  relax  machining  induced 
work  hardening  in  the  specimens.  Although  care  is  taken  t 
use  sharp  tools  and  slow  cutting  speeds,  no  annealing  is 
performed.  The  specimens  for  all  experiments  are  produced 
from  the  same  sample  of  aluminum;  and  the  uniaxial  stress 
strain  curve  for  each  specimen  is  assumed  to  be  the  same. 

From  the  least  complex  test  of  DPM-EPBEM ,  uniaxial 
tension,  one  expects  to  uncover  any  limitations  of  the 
hybrid  technique.  With  this  in  mind,  the  hybrid  analysis 
a  uniaxial  tension  specimen  is  described.  The  dimensions 
this  plane  stress  specimen  are  nine  inches  tall,  one-inch 
wide  and  .0625  inches  thick  (all  specimens  examined  have 
these  nominal  dimensions).  The  spots  are  applied  to  the 
specimen  as  indicated  in  Figure  19  and  as  shown  in 
Figure  20.  Two  interior  triangular  cells  are  used  to  mode 
the  region  enclosed  by  the  spots,  with  DPM  providing 
displacement  boundary  conditions  at  all  but  the  free 
surface.  One  should  note  that  only  the  region  encompassed 
by  DPM  is  modelled  by  the  boundary  element  method  and  that 
the  size  and  location  of  this  region  are  completely 
arbitrary.  Figure  12  shows  a  plot  of  the  x-component  of 
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Figure  19.  Schematic  of  the  spot  locatior.s  or.  the 
uniaxial  tension  specimen. 

displacement  which  is  calculated  'ey  OPM  at  (x=0.25  in.  , 
y=0.0  in.).  The  results  conform  to  the  shape  one  expects 
fer  a  load-displacement  curve  of  weakly  work  hardening 
materials.  A  little  more  convincing  is  the  plot  shown  m 
Figure  22,  which  provides  the  stress-strain  curve  calculated 
by  DPM-EPBEM  along  with  the  experimental  version.  Two 
interesting  characteristics  of  note  are  the  good  correlation 
between  the  MTS  and  DPM-EPBEM  results  and  the  limited  useful 
range  of  DPM-EPBEM. 
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initial  loading  steps  and  although  it  is  less  pronounced 
than  the  rigid  body  rotation  fr~m  specimen  misalignment, 
this  slippage  is  certainly  significant. 

DPM  is  able  to  discern  the  rigid  body  motion  from  both 
sources  and  passes  this  rigid  body  motion  information  to 
EPBEM.  Unfortunately,  EPBEM  expects  the  constrained 
displacements  and  has  no  means  to  separate  out  the  rigid 
body  motion.  Therefore,  the  boundary  element  routine 
calculates  erroneous  stress  and  strain  values. 

Whereas  DPM  is  responsible  for  the  lower  performance 
limit,  EPBEM  is  responsible  for  the  upper  performance  limit. 
EPBEM  expects  a  work  hardening  material  and  has  difficulty 
converging  to  a  solution  for  a  perfectly  plastic  material. 
The  upper  limit  results  from  the  fact  that  1100-H14  aluminum 
behaves  in  ai  almost  perfectly  plastic  manner.  At  high  a 
stress,  the  flat  portion  of  the  stress  strain  curve  causes 
the  iterative  solution  procedure  to  become  unstable; 
essentially,  the  initial  strain  procedure  fails  to  converge 
to  a  single  value.  It  is  believed  that  the  load  increments 
in  the  flat  portion  of  the  stress-strain  curve  produce 
strain  increments  too  large  for  the  initial  strain  solution 
procedure  to  assimilate.  An  upper  performance  limit  is 
typical  in  all  of  the  experiments  performed  but  it  varies 
with  loading  conditions  and  geometry.  Still,  the  DPM-EPBEM 
upper  performance  limit  usually  occurs  at  an  equivalent 
stress  of  16,500  psi.  Overall,  tne  uniaxial  tension  test 


30 


exhibits  dependable  results  within  a  specific  equivalent 
stress  range  and,  as  anticipated,  unnasks  basic  weaknesses 
of  the  DPM-EPBEM  hybrid  technique. 

IV. 5  Perforated  Strip  Tensile  Test 

The  uniaxial  tensile  test  produces  comforting  results. 
In  addition,  the  results  establish  both  the  accuracy  and  the 
performance  range  of  the  hybrid  technique.  This  next  test 
seeks  the  solution  to  a  more  realistic  engineering  problem 
and,  among  other  important  details,  the  perforated  strip 
exposes  DPM-EPBEM  to  a  stress  state  containing  high  stress 
gradients.  The  stresses  calculated  by  the  hybrid  technique 
are  compared  to  both  AN SYS  finite  element  and  Theocaris  and 
Marketos'  experimental  results  [95]. 

The  spots  are  applied  to  the  specimen  at  the  locations 
specified  in  Figure  23  and  shown  in  Figure  24,  where  the 
domain  enclosed  by  the  spots  is  moAoiied  as  seen  in  Figure 
25.  A  careful  comparison  between  Figure  25  and  Figure  23 
reveals  that  not  all  of  the  DPM  spots  and  EPBEM  nodes  are 
found  at  the  same  locations.  The  finite  size  of  the  spots 
establishes  a  practical  limit  to  the  spot  spacing; 
therefore,  DPM-EPBEM  uses  quadratic  interpolation  to  find 
the  displacements  at  nodes  which  do  not  coincide  with  spots. 


Figure  25.  Boundary  element  domain  discretization  of 
the  perforated  strip. 
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Figure  26.  Finite  element  grid  for  the  perforated 
strip . 

(x=0.25  in.,  y=0.0  in.)  as  a  function  of  load  step  are  given 
in  Figure  27.  From  here  on,  the  stress  concentration  factor 
(SC F)  is  defined  as  the  current  '/-component  of  stress 
divided  by  the  applied  stress.  The  finite  element  and  DPM- 
EPBEM  results  compare  veil,  with  ETM-EPBE.'l  exhibiting  one 
larger  excursion  at  the  sixth  load  step.  No  immediate 
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The  difference  in  the  response  could  be  caused  by  one 
of  many  reasons.  Theocaris  and  Marketos  use  a  bi-linear 
approximation  of  the  uniaxial  stress-strain  curve  in  their 
analysis  whereas  the  DPM-EP3EM  uses  a  combination  of 
Ramberg-Osgood  and  linear  fits.  Theocaris  and  Marketos  make 
no  mention  of  whether  they  did  or  did  not  annea"  their 
specimens,  nor  do  they  discuss  their  loading  apparatus  in 
any  detail.  Therefore,  one  would  not  necessarily  expect 
perfect  agreement. 

Other  interesting  information  produced  by  the  hybrid 
technique  is  the  load-displacement  curve  for  (x=0.Q  in., 
v=.25  in.)  which  is  given  in  Figure  23.  The  load- 
displacement  results  mimic  (in  shape)  the  load-maximum 
strain  curves  reported  by  Zienkiewicz  [97],  Although  not 
definitive,  they  certainly  are  supportive  of  DPM-EPBEM 
reliability.  In  his  study,  Zienkiewicz  used  the  material 
properties  reported  by  Theocaris  and  Marketos. 

A  final  plot.  Figure  29,  demonstrates  the  extent  of  the 
plastic  zone  by  giving  the  SC F  at  various  radial  locations 
from  the  root  of  the  hole  (for  ty=6000  psi) .  The  slight 
inflections  at  x=0.4  in.  and  x=0.275  in.  are  apparent  in 
all  of  the  graphs  provided  by  Theocaris  and  Marketos  (95' 
and  are  therefore  reassuring. 
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IV.  6  V-r.otched  Specimen  Tensile  Test: 

The  correlation  between  DPM-EPBEM,  finite  element 
analysis  and  the  purely  experimental  results  establish  a 
footing  from  which  a  more  complex  problem  can  be  approached. 
The  concluding  experiment  attempts  to  subject  the  hybrid 
technique  to  a  stress  state  truly  representative  cf  common 
engineering  problems,  a  notched  tensile  specimen  [98]. 

This  geometry  (Figure  30)  produces  high  stress  gradients 
localized  near  the  root  of  the  notch  with  a  complex  stress 
field  in  the  surrounding  material  [99].  It  should  be  noted 
that  there  is  a  small  machine-induced  rounding  of  the  notch 
root;  therefore,  the  notch  cannot  be  considered  sharp  in 
the  ideal  sense.  Once  again,  the  ANSI’S  finite  element  code 
and  selected  published  results  are  compared  to  those 
produced  by  the  DPM-EPBEM  hybrid  technique. 

As  before,  the  spots  are  positioned  according  to 
Figure  30  and  displayed  in  Figure  31.  A  sketch  of  the 
boundary  element  grid  can  be  found  in  Figure  32  and  one  cf 
the  finite  element  grid  in  Figure  33.  Again,  the 
homogeneous,  isotropic  hardening  option  and  two- 
dimensional,  four  point  isoparametric  elements  are  used  in 
the  finite  element  solution  of  this  problem. 

Both  the  SCF  vs  load  increment  (Figure  34)  and  the  load 
-displacement  curve  for  (x=0.0  in.  ,y=.47  in.)  (Figure  36) 
are  provided  for  evaluation  of  the  performance  of  DPM-EPBEM. 


The  finite  element  stress  concentration  factor  plct 
calculated  by  first,  appr/ir.q  an  end  load  of  2000  psi 
then  adding  load  increments  or  11  psi  until  reaching 
of  5400  psi.  In  Figure  34.  the  r.yfcrid  and  FEM  resul 
differ  markedly.  The  first  point  of  deviation  is  at 
yield.  DPM-EPBE.M  predicts  that  tr.e  first  yield  will 
at  an  applied  traction  of  I-F0C  psi,  whereas  ANSIS  pr 


Figure  31.  Photograph  of  the  spots  applied  to  the 
V'-notched  tensile  specimen 
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the  first  yield  at  3100  psi.  One  cannot  say  which  is  a 
better  solution  but  domain  discretization  does  have  an 
impact  on  the  comparison.  If  the  notch  is  perfectly  sharp, 
then  plastic  deformation  initiates  when  the  first  load  is 
applied.  The  deformation  field  is  so  localized  that  many 
internal  cells  placed  near  the  notch  tip  would  be  needed  to 
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Figure  33.  Finite  element  mesh  for  the  arid  for  the 
V-notched  tensile  specimen. 

discern  the  motion.  The  stress  calculated  at  each  cell  is 
an  average  representation.  Therefore,  at  the  notch  root, 
the  larger  cells  produce  smaller  average  values.  A 
comparison  shows  that  the  finite  element  grid  has  three  less 
cells  near  the  notch  root  than  does  the  boundary  element 
grid;  therefore,  it  is  reasonable  to  expect  that  the  first 
yield  predictions  by  BEM  and  FEM  will  be  different. 
Experimental  error  also  contributes  the  different  initial 
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Figure  34.  Plot  cf  stress  concentration  factor  versus 
load  step  for  the  V-notched  tensile 
specimen . 

yield  predictions  of  DPM-EPBEM  and  FEM.  The  erratic  hybrid 
SCF  curve  suggests  that  seme  of  the  DPM  displacements  may  be 
in  error.  If  this  is  so,  then  the  early  increments  vcul^ 
equally  be  in  error.  Regardless  of  the  reasons  for  the 
difference  in  the  predictions,  the  qualitative  similarity  is 
apparent  with  the  quantitative  results  differing  at  most  by 
14.9%. 

Other  important  considerations  deal  primarily  with  the 
difference  between  the  ideal  loading  conditions  in  the 
numerical  simulation  and  the  actual  loading  conditions  in 
the  hybrid  technique.  Examination  of  the  post  yield 
specimen  (Figure  35)  indicates  that  yield  did  not  occur 


Figure  ?5.  Phctoaraph  of  the  V-Notch^d  Tcncii® 

Specimen  Before  Loading  and  After  Failure 
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simultaneously  at  both  notch  roots.  This,  in  turn,  implies 
that  the  specimen  was  not  perfectly  aligned  in  the  MTS 
grips,  precipitating  a  non- symmetric  stress  state.  in  such 
a  case  as  this,  one  cannot  expect  the  FEM  solution  to 
duplicate  the  experimental  results. 

Although  no  comprehensive  experimental  investigation  of 
the  V-notched  test  specimen  is  available,  inferences  from 
results  presented  by  several  previous  experimental  and 
numerical  studies  lend  credence  to  the  results  presented 
here.  The  load  displacement  curve  given  in  Figure  36 
qualitatively  compares  well  with  a  similar  numerical  curve 
reported  by  Telles  [18]  (page  171)  and  with  the  experimental 
curve  reported  by  Findley  and  Drucker  [100]  (1/16  inch 

specimens).  In  both  cases,  the  material  is  an  aluminum 
alloy  which  is  different  from  the  material  used  here. 
Additionally ,  the  aspect  ratios  and  the  net  cross  sectional 
areas  of  their  specimens  differ  from  those  used  in  the 
present  experiments.  Taking  these  facts  into  consideration 
and  the  fact  that  EPBEM  reached  its  upper  performance  limit 
before  the  plastic  limit  load,  Figure  36  exhibits  the 
characteristic  linear  load-displacement  curve  in  the  pre¬ 
limit  load  regime.  Linear  behavior  from  a  plastically 
deforming  body  may  seem  contradictory,  but  in  this  case  it 
is  reasonable.  During  the  early  stages  of  plastic 
deformation  (as  is  the  case  with  DPM-EPBEM)  the  plastic  zone 
is  localized  in  the  region  near  the  root  of  the  notch.  The 
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Figure  36.  Load-displacement  plot  for  (x=0,y=.46)  cn 
the  V-notched  tensile  specimen. 
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*inn re  37.  Plot  of  stress  concentration  factor  versus 
x-position  on  the  V-notched  tensile 
specimen . 
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surrounding  elastic  material  confines  the  plastic  zone  and 
effectively  stems  any  large  plastic  deformation.  The 
regions  away  from  the  root  are  essentially  unaware  of  the 
plastic  flow  until  the  loads  are  high  enough  to  produce 
substantial  redistribution  of  the  internal  stress  state  of 
the  specimen  and,  in  turn,  a  larger  plastic  flow  [100]. 

Finally,  the  stress  concentration  factor  along  the  y=0 
axis  (Figure  37)  follows  the  expected  course  by  increasing 
as  one  moves  closer  to  the  root.  The  EPBEM  upper 
performance  limit  (5400  psi)  prevents  a  fully  developed 
plastic  zone  and  hence  a  flattening  out  of  the  SCF  curve  is 
not  present.  The  SCF  versus  position  plot  is  taken  for  an 
applied  load  of  4700  psi,  a  load  at  which  the  plastic  zone 
is  still  confined  to  the  notch  root. 
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SECTION  V 
CONCLUSIONS 


The  two-dimensional  numerical-experimental  hybrid 
technique  described  herein  combines  context-free  syntactic 
pattern  recognition  (DPM)  and  elastic-plastic  boundary 
element  methods  ( EPBEM)  to  produce  a  useful,  non-destructive 
stress  analysis  tool.  Within  a  certain  load  range,  all  of 
the  presented  DPM-EPBEM  results  compare  well  with  finite 
element  and  purely  experimental  (where  available)  solutions. 
The  DPM-EPBEM  solution  for  the  V-notched  tensile  specimen 
exhibits  a  more  erratic  response  than  the  responses  of  the 
other  two  experiments.  Since  the  V-notched  tensile  specimen 
represents  a  state  of  stress  that  is  as  complex  as  one  would 
expect  to  encounter  in  engineering  problems,  an 
understanding  of  the  results  of  that  test  is  paramount  to 
evaluating  DPM-EPBEM's  usefulness  as  a  stress  analysis  tool. 

DPM-EPBEM  works  well  in  load  ranges  where  rigid  body 
motion  is  negligible  when  compared  to  the  displacements 
produced  by  distortion.  Rigid  body  motion  in  the  initial 
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stages  of  loading  (approximately  90%  of  the  proportional 
limit)  degrades  the  accuracy  of  the  DPM-EPBEM  hybrid 
technique.  This  lower  accuracy  limit  can  be  improved  by 
incorporating  rigid  body  motion  compensation  in  DPM  or 
EPBEM. 

In  instances  when  DPM-EPBEM  returns  questionable 
results,  a  possible  question  might  be  raised  about  how  to 
guarantee  that  the  applied  spots  are  in  the  exact  locations 
as  specified  in  Figures  49,  53  and  60.  The  accuracy  of  the 
spot  application  technique  described  in  Section  IV  is 
limited  by  human  error.  In  most  cases,  the  accuracy  of  the 
spot  location  technique  is  approximately  t.001  inch. 
Surprising  as  it  may  seem,  mirroring  the  schematics  is  not 
that  critical.  The  boundary  element  code  does  receive 
displacement  information  from  DPM,  but  it  also  receives  and 
stores  the  locations  of  the  spot  centroids  in  the  undeformed 
image.  Using  this  information,  DPM-EPBEM  interpolates  (or 
extrapolates  as  the  case  may  be)  the  displacement 
information  to  the  BEM  node  locations.  This  is  not  to  say 
that  errors  are  not  introduced  by  misalignment  of  the  spots. 
If  the  spots  are  to  be  located  along  a  horizontal  line,  then 
misalignment  in  the  vertical  direction  can  play  a 
significant  role.  However,  DPM-EPBEM  is  very  forgiving  with 
respect  to  spot  placement  errors. 


Given  that  the  spots  are  perfectly  aligned,  the  ability 
of  DPM  to  measure  the  displacements  accurately  is  certainly 
paramount  to  the  success  of  the  hybrid  technique.  In  the 
images  used  in  this  investigation,  typical  values  for 
signal-to-noise  ratio,  gray  level  difference  and 
displacement  resolution  are  200,  100  pixels  and  0.0001  inch, 
respectively.  In  many  situations  there  are  regions  in  the 
tested  specimens  where  this  resolution  is  insufficient.  As 
an  example,  the  stress  concentration  factors  for  the 
perforated  strip  and  the  V-notched  tensile  specimens  are 
calculated  at  the  boundary  of  the  respective  cutout  roots. 
Therefore,  the  DPM  displacement  error  plays  an  important 
role  in  the  SCF  results.  Both  cases  include  localized 
plastic  zones  and  the  same  nominal  dimensions.  Because  of 
these  localized  regions,  one  anticipates  that  the 
displacement  in  the  region  near  the  cutout  roots  governs  the 
accuracy  of  the  DPM  displacements.  Minor  displacement 
errors  away  from  the  cutout  roots  do  not  have  the  influence 
of  displacement  errors  in  the  plastic  zone.  The  plots  in 
Figures  57  and  64  follow  the  above  reasoning;  if  the 
resolution  for  both  experiments  is  the  same,  then  only 
regions  of  small  displacements  (on  the  order  of  the 
resolution  limit)  can  account  for  the  difference  in  the 
behavior  of  the  SCF  plots.  This  conclusion  immediately 
draws  attention  to  the  displacement  in  the  regions 
surrounding  the  cutout  roots.  The  plastic  flow  in  the 
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perforated  strip  is  not  as  confined  as  it  is  in  the  V-notch 
specimen.  Therefore,  the  displacement  at  the  root  of  the 
circular  cutout  experiences  larger  and  more  measurable 
displacements.  The  displacements  near  the  V-notch  are  much 
more  confined  by  the  surrounding  elastic  region  which  holds 
the  displacements  nearer  to  the  resolution  limit  of  DPM .  As 
a  result,  the  SCF  plot  for  the  V-notch  specimen  is  much  more 
erratic . 

One  should  note  that  the  specimens  used  in  the 
evaluation  of  DPM-FPBEM  are  smaller  (nominal  dimensions  are 
one  inch  wide,  nine  inches  tall  and  .0625  inch  thick)  than 
one  would  expect  to  encounter  in  many  true  engineering 
problems.  The  displacements  in  larger  specimens  will  be 
bigger  and  more  measurable  than  the  displacements 
encountered  here.  This  argument  holds  for  localized  plastic 
zones;  they  too  will  encompass  larger  areas,  providing  more 
sizable  displacements.  Therefore,  in  larger  structures,  the 
DPM-EPBEM  is  less  likely  to  encounter  the  resolution  limit 
of  DPM.  Avenues  for  improving  the  DPM  resolution  are 
available  for  problems  dictating  measurements  on  the  order 
of  the  resolution  limit.  The  interpolation  of  the 
displacement  information  from  spot  locations  to  node 
locations  has  the  effect  of  smoothing  out  the  influence  of 
the  poor  data  near  the  roots  of  the  cutouts.  This  suggest 
that  smoothing  techniques  tailored  to  these  types  of 
problems  may  be  useful.  The  real  solution,  however,  is  to 
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improve  the  resolution  of  DPM.  This  can  be  accomplished  by 
increasing  the  calibration  constants  through  magnificat:. dr. 
or  with  better  hardware,  or  by  compensation  for  system 
errors . 

Another  area  which  needs  improvement  deals  with  the 
elastic-plastic  boundary  element  method.  The  convergence  c 
the  initial  strain  elastic-plastic  boundary  element  sol-tic 
is  sensitive  to  load  increment  step  size.  As  seen  in  the 
experiments  presented  in  Section  IV,  it  is  not  always 
possible  to  apply  the  loads  in  such  way  as  to  stay  within 
the  realm  of  incremental  plasticity  and  thereby  guarantee 
convergence.  Therefore  a  solution  procedure  which  is  less 
sensitive  to  increment  size  would  be  more  desirable. 
Fortunately,  the  initial  stress  solution  approach  has  -use 
this  quality  [101].  It  is  anticipated  that  an  initial 
stress  solution  procedure  would  greatly  enhance  the  load 
range  in  which  this  hybrid  technique  is  applicable. 

The  evaluation  of  DPM-EPBEM  is  given  in  the  above 
paragraphs  so  that  the  limitations  of  this  technique  can  be 
understood  and  so  areas  of  further  research  can  be 
pinpointed.  But  even  in  its  present  stage  of  development, 
the  hybrid  experimental-numerical  procedure  described  in 
this  report  is  a  useful  stress  analysis  tool.  The  good 
results  presented  in  Section  IV  are  one  indication  of  the 
promise  that  the  DPM-EPBEM  procedure  holds.  A  second 
indication  of  the  promise  of  this  technique  is  that  it  is 
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fully  automated.  Finally,  the  equipment  required  to  produce 
a  hybrid  method  like  the  one  described  is  standard  at  most 
engineering  and  research  facilities.  By  incorporating  the 
above  suggestions,  the  displacement  pattern  matching, 
elastic-plastic  boundary  element  hybrid  method  can  be  shaped 
|  to  provide  good  solutions  over  an  even  wider  load  range. 
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APPENDIX 


The  two  volume  integrals  encountered  in  elasto-plast ic 
boundary  element  methods  appear  in  the  boundary  constraint 
equation  (Equation  11-15)  and  the  interior  strain  equation 
(Equation  11-40) .  Here  these  integrals  are  identified  and 
their  semi-numerical  integration  is  described. 

The  constraint  equation  volume  integral  is  written  as 

A2ijk  ~  sijk(r)dv 

J  v 

where 

sijk(r)  =  [C2(5ijr,k  +  5ikr,  j  "  <5jkr,i)  *  r,  i^ ,  j  r ,  k]  ( A-l ) 

with 

Cl  =-  ~ -  r 

4»  (l-i/) 


C2  =  (W)  , 
r  j_  =  -  f,  i  and 
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A  transformation  to  cylindrical  coordinates  simplifies 
the  integration. 


r*b  fR^>  * 

A2ijk  (r)  =  lira  1  1  sijk  (r,  j>  >  rdrdo  (A-2) 

J.pa  Je 

where  r,  <s>  .  a ,  R(j>)  and  lijk  are  defined  in  Figure  33. 

The  integration  is  carried  out  with  respect  the  line 
segment,  D,  perpendicular  to  R(c).  This  introduces  the  new 
coordinate  system  r1/  By  defining  the  derivatives  of  r 

in  cylindrical  coordinates,  Tijk  is  properly  recast. 


r 


-  j 


3  r  _  3  r  ■  ~  iL  +  3r  3  c 
3  Xi  3  r  i  3  X  a  3  ,'2  3  x  ^ 


( A-  3  ) 


As  seen  from  Figure  33,  =  cosj  ,  —  =  sino,  " — 1  =  e]_j_ 

3  T  i  d;  2  3  x  j_ 

and  ^ “ 2  =  e2i-‘  where  ej_j  are  the  direction  cosines  between 
3  x  j_ 

the  x  and  f  coordinate  systems.  Therefore, 

r ( j_  =  cos 4>  e]_i  +  sin<j>  e2i,  (A-4) 

and  combining  Equations  (A-l),  (A-2)  and  (A-4)  gives 
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Pb 

+  2  Ci  i  [  R  p  )  ( coso  ©ix  +  sini>e2i)  (cos^e^j  +  sinoe2j) 

J  a 

+  (cos,jexk  +  sini  e2k)]d$. 


From  the  geometry  shown  in  Figure  28,  j  =  « +<* ,  dj>=d;! 
and  Rp)  =  D/cos<a  ,  and  therefore  the  integral,  <tj_jk  is 


*ijk  =  Dclc2  1  (elk5  i j  +  elj,sik  "  eli5jk) 

+  tan«(e2k6ij  +  e2jSik  -  e2i6jk)) 

+  DC2  {  2cos25  e2_ie]_  je]_k 

+  2cososin<i  (e2ie1jelk  +  eiie2jelk  +  e1ie1je2k) 
+2sin2? (e2ieije2k  +  me2je2k  +  e2ie2jelk)  (A-6 

+  2sin2*tantfe2ie2je2k} . 


Upon  integration  Equation  (A-5) 
over  an  arbitrary  triangle. 


reduces  the  integral  of 


ASijk  =  C1C2D<  (elk<5  ij  +  eij<$ik  -  e1j_5jk)I1 

+ (e2k5 i j  +  e2 js ik  “  e2i«jk)I2)  +  DC2eliel jelkI3 
+  DC2 (e2ieijelk  +  j  e2n  elk  +  elieije2k)I4 

+  DC2(e2ieije2k  +  elie2je2k  +  e2ie2jelk)I5 
+  DC2e2ie2je2kI6, 


(A-7) 


where 


II  = 


*b 


d6  =  a  ' 


‘e  b 


12  =  I  tantfdtf  =  -In (COS-b)  =  ln(  b)  , 
Jfl,  cost)  a  r2 


r^b 


13  = 


J  a 


2cos^$d 6  -  (0b~^a)  +  (sin20  b~si*n2£  a)  ' 

2 


14  = 


^  b 


2cos0sin£d(J  =  sin2tf  b-sin2S  a  / 


f*b 


15  = 


2sin26d£  =  (Sb-Sa)  +  i(sin29b-sin2S a) 

2 


r#b 


16  = 


2tan0  Sin20  d£  =  21n(Eb)  +  cos2»b-cos29( 


The  volume  integral  associated  with  the  total  strain 
equation  is  simply  the  derivative  of  Equation  (A-7)  (see  the 
development  of  Equation  11-40),  that  is 

AE^ijfc  =  —  (AZijfc) .  (A-8) 

dZ  2 

Here  the  derivative  is  taken  with  respect  to  the  load  point 
If  f  is  perturbed  then  D  and  its  associated  angles  change; 
C]_,  C2  and  e^j  remain  fixed.  Evaluating  Equation  (A-8)  is 
now  reduced  to  finding 


'T"  " 


109 


Ij^  =  - - (DIi*  where  l  =  1,  2,  ...  6, 

3<  2 


which  by  the  chain  rule  expands  to 


Iii  = 


3D  Ij_  dj_ 

3 <  2  3  8  3^2 


(  a-9  ; 


3D  ^9 

The  first  term  - —  =  -eii  and  in  the  second  term  is 

3  £  ^  3  £  i 

evaluated  as  follows: 

8  =  tan-1  (^2)  ~  a  , 

•  1 


<LL  _  _  Lz 

3  (  R 


( A- 10 ) 


but  ;■  2  =  Rsin(S+<x),  after  expanding  and  substituting  into 
Equation  (A-10) 


d  8  _  sin<cosa  -  costs 

at  ,  "  R 


8  sma 


similarly 


2X  cosjcosa  -  sinasim 

3*,  '  R 


By  making  use  of  the  direction  cosines 


U.  -  .  e2iC°S6 


-  &l£Sln9 


R 


110 


Now,  by  substituting  the  above  results  into  Equation  (A-9) 
AElijk  can  be  written  for  a  representative  triangle 
AEiijk  =  clc2(elk*ij  +  el  j15  ik  "  eii^jk)1! 

+  C1C2(e2k6ij  +  e2j*ik  “  e2i5jk)I2 

+  C2elie1jelkdl3  +  C2(e2ie1jelk 

+  e1ie2jelk  -  eiiSlj^k)1-; 

+  c2 ( e2 iel j  e2k  +  elie2 j  e2k  +  e2ie2jelk)x5 
+  c2e2ie2je2k  1 6  - 

where 

X1  =  ~ ( s b_tf  a) ell  +  (ib'ia)  * 


I2  =  v btan!?b  "  iatan!,a  ~  ln(~b)e12  , 


13  =  2rjaCOS29a  -  2r>bCOS2$b 

-  {(«b_t,2)  +  sini?bsinf?a  -  cosO  asin9  a)  e12 

14  =  2r7bsini?  bcostf  b  -  2n  asini?  acosS  a 

-  (sin20b  -  sin29a)e12 

15  =  2r?bsin2Sb  -  2  r;  asin2^  a  -  (^b_sa)e12 

+  (sin0bcosflb  -  sintf acos0 a) e12 


2r?btan5  bsin20  b  -  2r,atan0  asin20  a  -  21n(£a)e12 

,  2  2  ra 

-  (cos^b  -  cos^0a)e12; 


d  S 

where  r?a  and  r;b  are  D  evaluated  at  points  a  and  b, 
respectively. 
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